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DFT Discrete Fourier transform 
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1 Introduction 


1.1 Ultrasonic time delay measurement systems 


Most ultrasonic measurement technologies rely on the time delay esti- 
mation (TDE) of ultrasonic pulses. Depending on the context, the time 
delay is also known as time-of-flight or transit-time. In other words, an 
ultrasonic transducer generates a pulse, which propagates through a 
medium and gets either reflected back to the transducer operating in 
sensing mode or received by a different transducer. The resulting time 
delay can then be estimated and used to infer the desired measurand, 
e.g., the speed of sound (SoS), the distance to an object, or the veloc- 
ity of a medium. Typical applications include, but are not limited to, 
non-destructive testing (NDT) [66], sonar [1], parking sensors [125] and 
ultrasonic flow measurements (UFM) [79]. As an example, figures 1.1(a) 
and 1.1(b) show the different principles of parking sensors and UFM, re- 
spectively. While UFM requires two transducers on opposite sides of the 
pipe, parking sensors can be realized using a single transducer operating 
alternatively as a transmitter and receiver. 


1m a pee 


w, Re 
mr S77 
(a) Ultrasonic flow measurement. (b) i ee sensors. 


Figure 1.1 Ultrasonic transit-time estimation applications. 


As there are a myriad number of possible ultrasonic applications, the 
application-specific challenges are equally varied. However, all systems 
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have in common that their accuracy is mainly determined by the mea- 
surement of the underlying time delay. The challenge of an accurate 
TDE can be approached in two ways. Either the signals are improved 
by adapting the design of the ultrasonic system, or the signals are sub- 
sequently enhanced by using advanced signal processing methods—a 
combination of both is also possible. Firstly, the possibilities by adapting 
the design are considered. As the design of an ultrasonic measurement 
system contains a lot of degrees of freedom, the following list outlines a 
few key points that can be varied: 


= The measurement principle represents the most important aspect 
in the design to get a strong measurement effect. While NDT for 
hidden defect detection, sonar, and parking sensors are all based 
on the pulse-echo method [1, 8, 17], other principles are using the 
resonance of a system for TDE or the superposition of the SoS 
with the velocity of the medium. For example, Figure 1.1(b) shows 
that the ultrasonic waves transmitted by the sending transducers 
are reflected at other cars. The echo can be recorded either by the 
sending transducer, which can be switched to operate alternately in 
sending and receiving mode, or by another receiving transducer. In 
contrast, as shown in Fig. 1.1(a), UFM is based on the superposition 
of the SoS with the velocity of flow (VoF). Here, waves sent in 
upstream direction are slower than in downstream direction. 


= The operating frequency range is chosen depending on the medi- 
um and the precision requirements. Higher frequencies are tech- 
nically more expensive and lead to improved accuracy, but due 
to the larger attenuation at higher frequencies, the signal-to-noise 
ratio (SNR) deteriorates. If the medium is air, the used frequency is 
typically limited to 50 kHz, because the attenuation at high frequen- 
cies is significantly larger than in liquid media. Typical frequency 
ranges for the diverse applications are summarized in Table. 1.1. 


= Finally, the choice of transducer type depends on the desired fre- 
quency range and the type of ultrasonic waves to be generated. 
Nowadays, most transducers are based on piezoelectric ceramics, 
which convert an electrical charge to a mechanical displacement 
and vice versa. They can be arranged in arrays and combined with 


1.2 Influence of interfering signals 


Table 1.1 Typical frequency ranges for a selection of ultrasonic applications. 


Application Frequency range 
Parking sensors 40 — 50 kHz [91] 
Maritime sonar (mid-frequency) 2.6 — 50 kHz [22] 
NDT (air-coupled) 0.75 — 2 MHz [8] 
UFM (liquids) 0.5 — 4 MHz [103] 
UFM (gas) 200 — 500 kHz [11] 
Medical ultrasound imaging 2 — 15 MHz [51] 


beamforming to reach better directional characteristics [1], or sim- 
ply be connected with backing and electrodes and then placed 
inside a metal case. In the case of non-invasive transit-time UFM, 
so-called clamp-on sensors [114] that generate surface acoustic 
waves are often preferred. 


While the development and improvement of ultrasonic systems re- 
quire a lot of iterations and are, therefore, very time-consuming and 
expensive, the digital signal processing method can be easily adapted 
using a software update, even in the late stage of the development. Hence, 
the focus of this thesis lies on signal processing techniques for TDE in 
ultrasonic systems. Furthermore, this thesis places a special focus on the 
TDE accuracy and robustness against interfering signals and noise, as 
these properties are usually of high importance. 

From a signal processing point of view, ultrasonic time delay measure- 
ments can be simplified into the question, how to get an accurate TDE 
from a 1D measurement signal sq (t) and a reference signal s,(t). Besides 
signal distortion due to the acoustic transmission system, the most likely 
reasons that complicate the TDE are unwanted signals, such as noise or 
interfering signals, which are discussed in more detail in the next section. 


1.2 Influence of interfering signals 


Unwanted signals in ultrasonic measurements can be separated into 
external disturbances due to uncorrelated sound sources, interfering 
signals due to multipath propagation, and stochastic noise such as addi- 
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Sending wedge transducer Receiving wedge transducer 


Pr Rayleigh wave 
IA 


Breaking crack 


[Testing material 


Backwall reflection 


Figure 1.2 Interfering signals (red) due to multipath propagation in an NDT scenario. 


tive white Gaussian noise (AWGN). A typical scenario with unwanted 
signals due to multipath propagation is depicted in Figure 1.2, where a 
homogeneous material is tested for cracks or other defects. Due to the 
impedance jump at the crack, the incident wave gets reflected towards 
the receiving transducer, which in turn allows locating the crack by esti- 
mating the time delay, if the SoS of the material is known. However, the 
signal created by the crack, namely the target signal, is superposed by a 
direct propagation between the transducers along the surface of the test 
specimen and by back wall reflections. An external disturbance could 
be caused by sound sources in the test environment, which makes them 
uncorrelated to the target signal. Additionally, as with any measurement, 
the recorded signals are subject to stochastic noise, caused by polariza- 
tion noise or thermal noise of the piezoelectric ceramic as well as the 
thermal noise of the measurement electronics. 

While the influence of AWGN can be reduced by averaging and uncor- 
related sound sources can be suppressed by using modulated excitation 
signals, the separation of interfering signals is more difficult. Since they 
are highly correlated with the target signal, a modulated excitation signal 
shows the same effects on both the interfering and target signals, and 
therefore introduces no information to separate them. Furthermore, the 
interfering signals are in the same frequency range as the target signal 
and stationary in the sense that they do not vary over multiple measure- 
ments, which renders frequency filtering or averaging useless. If a short 
pulse is used as excitation and the transducer bandwidth is sufficiently 
high, the different signals can be separated by gating in the time domain 
(see Fig. 1.3(a)). Otherwise, as it can be seen in Fig. 1.3(b), the interfering 
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|. Interfering Signals | H Interfering range =] 
Am AN An 
l l | | | | 
40 50 60 70 40 50 60 70 
tin ps tin ps 
(a) Separable interfering signals. (b) Non-separable interfering signals. 


Figure 1.3 Superposition scenarios of interfering signals (the measurement signals are 
colored in blue and the target signals in red). 


signals are located in the same time and frequency range resulting in a 
distorted observed signal, in the literature often termed as non-separable 
interfering signals [69, 102]. 

The main problem arising from the interfering signals is their influence 
on the accuracy of the TDE. Since they show no measurement effect, the 
actual time delay of the target signal may be hidden under the time delay 
of the interfering signals, which induces a systematic estimation error. 
Depending on the estimation method, the error can be significant, if 
for example the target signal cannot be detected at all and the peak of 
the interfering signals is taken for the time delay of the target signal. 
Even if the level of the interfering signals is small compared to the target 
signal, the accuracy is considerably reduced. To illustrate this, a thought 
experiment of a TDE in the presence of interfering signals is conducted 
in the following. For simplicity, pure sine signals are used as the target 
and interfering signals, which is at the same time the worst case, because 
the bandwidth of such a measurement is nearly zero. The time delay gets 
estimated via the phase shift of the measurement signal since it describes 
the time delay for pure sinusoidal signals. To describe the problem, the 
measurement signal 


salt) = s(t — 7) + 5; (t) (1.1) 
is composed of the target signal s,(t) and additive interfering signals 
s,(t), which are defined by 

s(t) = A, sin(wot), (1.2) 

s;(t) = A; sin(wot + p). (1.3) 
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I 
1/fy -| ——ScNR = 10dB 


— ScNR = 5dB 
— ScNR = 3dB 


—— ScNR = 1dB 


0 r 1/fo 


Figure 1.4 Worst-case estimation error caused by the interfering signals with y = 1/2 
depending on their level (ScNR = 20 1g(A,/A,)dB). 


Given the equations (1.1)-(1.3), the observable time delay of the measure- 
ment signal can be calculated using the harmonic addition of sine waves 
for arbitrary phase shifts [12, p. 84] resulting in the estimated time delay 


( A, sin(woT) + A; sin(p) ) 
A, cos(wyr) + A; cos(p) 


(1.4) 


2 1 
T = — arctan 
Wo 


The results of this equation can be seen in Figure 1.4, where the esti- 
mated time delays for different signal-to-correlated-noise ratios (ScNR) 
are shown. This ratio must not be confused with SNR, as ScNR describes 
the level of interfering signals and not the stochastic noise. It can be seen 
that the error, oscillating with fo, is dependent on the real time delay 7 
and the ScNR. The absolute maximal deviation 7 —7 is determined by the 
ScNR, while the phasing is determined by y, which has been set to 7/2 
in this example. Even for small levels of interfering signals, where the 
amplitude of the target signal A, is about three times higher than the in- 
terfering signals amplitude A; (ScNR = 10dB), the estimation error can 
be up to 5% of the period duration 1/ fọ. In summary, the errors induced 
by the interfering signals are systematic, dependent on the phasing of 
the superposition and can only be mitigated by increasing the excitation 
frequency or damping the interfering signals mechanically. 

Another source of error that occurs in the presence of interfering sig- 
nals is the influence of process parameters, e.g., temperature or pressure. 
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Since the estimation error is dependent on the superposition, any process 
parameter that leads to a deviation of the SoS changes the phasing and, 
therefore, the estimation quality. 

One scenario, where interfering signals are especially problematic, 
is ultrasonic flow metering, e.g., using clamp-on sensors, as shown in 
Fig. 1.1(a). Here, the ultrasonic waves are transmitted into the fluid via the 
pipe wall, whereby the interfering signals account for a significant pro- 
portion of the measurement signals. The different approaches reported 
by the literature to address this problem are based on either mechanical 
damping or using flow statistics. However, their constraints limit their 
applicability, since mechanical damping is usually temperature- and 
frequency-dependent, and the assumed flow statistics are too restrictive. 
Therefore, extended algorithmic methods for filtering the interfering 
signals are proposed in this thesis, exploiting the property that the target 
signals exhibit different dynamic behavior compared to the interfering 
signals over multiple measurements. In contrast to the state-of-the-art 
(SotA) methods, the presented algorithms require little dynamic behavior 
of the target signals and also work over a wide temperature range. 


1.3 Scope of the thesis 


This thesis presents model-based algorithmic approaches for interference- 
invariant TDE, which are specifically suited for the estimation of small 
time-delay differences with a necessary resolution well below the sam- 
pling time. Therefore, the methods can be applied particularly well for 
transit-time UFM, since the problem of interfering signals is especially 
prominent in this application. The main focus lies on how to use multiple 
measurements with varying time delays or process parameters for the 
separation of interfering signals in UFM, while maintaining good robust- 
ness against AWGN and a high resolution. To this end, a signal model 
is assumed that contains stationary interfering signals, which are not 
dependent on changing time delays, and a target signal, which contains 
the measuring effect. 

Firstly, the signal model of an UFM and its dynamic behavior under 
temperature or time delay variations is examined in Chapter 3. The aim is 
to generate valid simulated data sets, which are used to test the developed 


1 Introduction 


methods both under the premise that the data fits to the signal model 
perfectly and under the premise that model errors are present. In the 
course of this, the properties of the signal model components, such as 
bandwidth, stationarity and temperature dependency, are identified. For 
this purpose, a new method to model the temperature dependence of 
the interfering signals is presented. After the characterization of the total 
measurement system, the signal model—adapted to UFM—is used as 
the basis for two new methods whose objective is to reduce the impact 
of the interfering signals. 

The first proposed technique, described in Chapter 4, extends the dy- 
namics-based approaches of the literature by weakening the constraints 
on the needed variance of the time delays. To this end, a new representa- 
tion of multiple measurement signals as point clouds is introduced. The 
point clouds are then processed using the principal component analy- 
sis (PCA) and B-splines, leading to either interference-invariant TDE or 
estimated interfering signals. In this context, a novel joint B-spline and 
misalignment approximation is developed to enhance the robustness. 

The second approach consists of a regression-based estimation of the 
time-delay differences by learning reasonable signal subspaces. These 
subspaces are efficiently calculated by the analytic wavelet packet trans- 
form (AWPT) before the resulting coefficients are transformed into fea- 
tures that correlate well with time-delay differences. Furthermore, a novel 
subspace training approach that works unsupervised is proposed and 
compared to the conventional filter- and wrapper-based feature selection 
methods. 

Finally, both methods are tested in an experimental UFM system with 
a high level of interfering signals present, where it is shown that they 
are in most cases superior to the literature methods. The quality of the 
methods is evaluated using the accuracy of the TDE, since the ground 
truth for the interfering signals cannot be determined reliably. Different 
data sets are used to analyze the dependencies on the hyperparameters, 
the process conditions and, in the case of the regression-based method, 
the training set. 


2 State of the Art 


In this chapter, an overview of the relevant literature concerning TDE, 
transit-time UFM, and suppression of interfering signals is given. The 
first section presents and categorizes SotA methods for TDE that will be 
used in Chapter 6 as a reference for the newly developed methods. While 
many of the methods are only described in a brief overview, the methods 
that are to be used as a reference are explained in more detail with the 
necessary equations. Since the focus of this thesis is the robustness against 
interfering signals using the transit-time UFM as an application example, 
sections 2.2 and 2.3 outline the SotA regarding transit-time UFM and 
interfering signals suppression, respectively. 


2.1 Time delay estimation 


In TDE problems, the objective is to find the time delay between a refer- 
ence signal and its delayed and attenuated echo in the presence of noise. 
Let the general problem be modeled by 


s,(t) = s(t) +n, (t), 


(2.1) 
salt) = A, s(t — 7) + na(t), 


with the reference signal s,(t), the delayed signal s4(t), the target signal 
s,(t), the attenuation A,, and the time delay 7 [61]. In most cases, the 
AWGN signals n,(t) and n4(t) are assumed to be uncorrelated with each 
other and with s,(t). 

The different estimation methods that are reported in the literature can 
be categorized into three classes according to Fang etal. [30]: time-based 
[17], correlation-like [5] and waveform-fitting [3, 47, 87]. While time- 
based methods are easier to implement and require little computational 
effort, correlation-like methods yield better accuracy and robustness 
against noise for low SNR [117]. If prior knowledge such as propagation 
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Figure 2.1 Zero crossing estimation via linear approximation of the Hilbert angle &,. The 
zero crossing to be selected must be defined by another criterion. 


characteristics or transducer impulse response is available, the waveform 
can be physically or empirically modeled to estimate the time delay 
via an optimization approach. However, the necessary parameters are 
often only known approximately which reduces the accuracy. For this 
reason and the significant computational effort of physically modeling 
the acoustic waves, waveform-fitting approaches are not widespread. 


2.1.1 Time-based estimation methods 


Let us first consider time-based methods, which are based on the determi- 
nation of prominent points in the measurement signal, e.g., intersections 
with thresholds [35], zero crossings [139] or envelope edges [4]. 

As early as 1974, Frederiksen and Howard [35] implemented a thresh- 
old-based determination of the time of arrival of wave packets on a 
monolithic chip. Via a specially developed analog circuit, an impulse 
noise rejection has even been integrated, which discriminates against re- 
ceived signals that are too short. Since a fixed threshold is not very robust 
to fluctuating signal intensities, various extensions of the threshold-based 
TDE to multiple thresholds [71] or variable thresholds [140] have been 
introduced. However, thresholds cannot only be applied directly in the 
time domain. For example, Duarte et al. [28] have presented the use of 
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thresholds in spectrograms at a given frequency as a way to perform 
TDE. 

Alternatively to the threshold-based TDE, the zero crossing can be used 
to determine the time delay—an easy algorithm for the zero-crossing 
detection is given by Zhou etal. [139]. To get a more robust estimation of 
the zero crossing, Kupnik etal. [63] use the Hilbert transform in combi- 
nation with a linear approximation of the phase signal. The approach 
is depicted in Figure 2.1. As it can easily be seen, narrow-band signals 
lead to a linear behavior of the phase signal ¢,,(¢), which results from 
calculating the argument of the analytic signal s,,(t). The robust zero 
crossing is then extracted as the intersection of the linear phase fit Da (t) 
with the time axis. The same approach is also presented by Roosnek 
[99], where the linear approximation ¢,,(t) is calculated by the weighted 
least squares using the signal envelope |s 4 as(t)| as weights. This specific 
method—in the following simply called the zero-crossing method—is 
one of the SotA methods that are implemented as a reference to measure 
the accuracy of the methods proposed in this dissertation. It has to be 
noted that the zero-crossing method is ambiguous since in a sinusoidal 
signal there are usually several zero crossings. A solution to always de- 
termine the time delay of the same zero crossing is proposed by Fang 
etal. [30], where finding the same zero crossing is guaranteed by first 
estimating the signal onset using a classification. 

In addition to the methods already listed, many more elaborate meth- 
ods exist in the literature. They are based on, for example, features such 
as the intersection of an envelope tangent with the time axis [4], or they 
determine the time delay by tracking multiple sample points individually 
[136]. 


2.1.2 Correlation-like time delay estimation 


Correlation-based methods for TDE have been published since the early 
1970s, as TDE is one of the main applications for cross-correlation. The 
core of each method lies in the cross-correlation function 


Ra, = E{s,(t) salt + 7)} 


x a s,(t) sat + 7)dt, a 


co 
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which shows a peak at the real time delay 7 under ideal conditions [61, 89, 
93]. A simple application of this function is the cross-correlation flowme- 
ter, where two axially separated transmission paths are used to measure 
the movement of bubbles, turbulence, or other clumps of particles [5]. 
However, the ideal conditions include the noise signals being indepen- 
dent zero-mean stationary Gaussian and the available signal sequences 
being ergodic and long enough to approximate the expectation value by 
time integration. In real applications, the observation time, the sampling 
rate and the bandwidth of the target signals s,(t) are usually limited, 
leading to a noise polluted spread out cross-correlation. Moreover, mul- 
tipath propagations are neglected in the standard approach which can 
further deteriorate the estimation quality. To cope with these problems, 
several extensions of the classical method have been developed and 
published over the years. These extensions range from pre-filtering the 
signals to advanced approximation techniques based on the calculation 
of weighted cross-correlations in the Fourier domain. 

An important class of cross-correlation techniques that needs to be con- 
sidered is the generalized cross-correlation. This class is characterized by 
the fact that the cross-correlation calculation is carried out in the Fourier 
domain with a frequency weighted spectrum. Knapp and Carter [61] 
compared different frequency windows H, |, (f) for the generalized 
cross-correlation method, in which the time delay i is estimated by 


= arg max i. (7) 


=argmax | Hy, a (f) Sa(f) SA) af, (2.3) 


where S,(f) and S,(f) denote the Fourier transforms of s,(t) and s,(t), 
respectively. The aim is to improve the estimation quality of the cross- 
correlation by emphasizing frequencies with high SNR. This approach is 
equivalent to pre-filtering the signals with h, (t) and h, (t), respectively, 
before the cross-correlation is calculated. Since the choice of H, |, (f) 
influences both the robustness and accuracy of the estimation, a selection 
of different weighting schemes, also known as processors, is listed below: 


= Roth’s processor 


E, 


S.A) SA (2.4) 
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suppresses frequencies with low SNR but also spreads out the 
cross-correlation en (r) [101]. 


Using the phase transform processor 


u 1 
BRETT 


the cross power spectral density is normalized, which can sharpen 
the peak of the generalized cross-correlation [16]. 


A, sP) 2.5) 


The Eckart filter 


SA SA) 
NAP) Nef): NaP NIO) 


has been developed by Eckart [29] as a filter that maximizes the 
deflection. In this context, deflection represents the ratio between 
the filter’s response to the target signal compared to the response 
to noise. It also has the property to suppress frequency bands with 
low SNR, but the spectral densities of the target signal and the 
noise have to be known beforehand. 


A, (f= (2.6) 


The last spectral filter to be mentioned is the Hannan-Thomson 
processor [42], for which the approximation 


u ISA 
INCA? INE AP INA Nah)? 
was proposed by Brandstein et al. [10]. Note that this filter leads to 


the maximum likelihood estimate of the time delay, with an error 
variance approaching the Cramér-Rao lower bound (CRLB) [61]. 


Hy Ad) (2.7) 


All above mentioned generalized cross-correlation methods need the 
power spectral densities (PSD) of either the target signal, the noise signals, 
the measurement signals, or all together. In practice, these densities are 
approximated via the periodogram, Bartlett’s/Welch’s method, or other 
spectral density estimators. Alternatively, a data-driven FIR filter design 
approach, which does not rely on spectral densities, is presented in a 
newer work [6]. Nevertheless, as already stated above, all methods are still 
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limited by the bandwidth and the observation time, leading to a smeared 
cross-correlation function with additive residual noise. If the peak of the 
cross-correlation function is too broad or polluted by noise, Cabot [13] 
has shown that the Hilbert transform can then be used to convert the 
maximum of R, ,, (7) into a zero crossing. This technique could also be 
combined with the pre-filters in the generalized cross-correlation. 

In contrast to using the peak of the cross-correlation, other proper- 
ties such as the phase of the cross power spectrum can also be used 
to estimate the time delay, especially in cases where the noise signals 
na(t) and n,(t) are correlated with each other [93]. Since correlated noise 
sources are a typical scenario and originate from external sources in the 
environment, this problem was also addressed by Nikias and Pan [89] 
by employing higher-order spectrum domains such as the bispectrum. 
Their approach is based on the property that the third moment of a 
zero-mean Gaussian process vanishes, which removes the influence of 
the correlation between the noise signals. However, the term correlated 
noise must not be confused with the interfering signals, the reduction of 
which is the objective of this thesis, as these are not only correlated with 
each other but also with the target signal. 

While all methods based on the cross-correlation perform TDE by max- 
imizing an adapted form of (2.2), an alternative correlation-like technique 
is based on the minimization of the average square difference function 
(ASDF) 


RASPF (7) = i P (s(t) — salt + N)? (2.8) 


co 


or the average magnitude difference function (AMDF) 


RISMEG) = | leE) satt + Pleat, (2.9) 
which were examined in the works of Jacovitti and Scarano [50] and Fert- 
ner and Sjolund [31]. These techniques are also known in the literature as 
the sum-of-squared-differences and sum-of-absolute-differences method, 
respectively [36, 121]. Compared to the conventional cross-correlation, 
ASDF and AMDF require less computational effort since there is no need 
for multiplications and normalization. Furthermore, the variance of the 
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estimation error is closer to the CRLB than the cross-correlation in a 
white signal scenario [50]. 

A special challenge that has been neglected so far in the cross-corre- 
lation methods is multipath propagation. If the difference of the time 
delays between the system’s various propagation paths is large enough, 
the cross-correlation shows distinct peaks that correspond to the time 
delays of the different paths. The limiting factor for the necessary differ- 
ence between the time delays is the available signal bandwidth which 
is inversely proportional to the width of the correlation peaks [55]. For 
high SNR, the resolution can be improved by whitening the spectrum 
as shown by Khyam etal. [55]. The other approaches consist of using 
pulsed-wave ultrasound as excitation in combination with a windowed 
cross-correlation, where the signals are windowed in the time domain 
before calculating the cross-correlation [9, 129]. While fixed rectangular 
windows are often used [33, 45], Xu et al. [133] presented an algorithm 
that calculates the cross-correlation in the wavelet domain with care- 
fully selected scaling levels. A closer look into this algorithm leads to 
the conclusion that this is a special case of a time-frequency windowed 
cross-correlation showing similarities to a time-windowed generalized 
cross-correlation. 

The last step of each correlation-like approach is the subsample approx- 
imation of the time delay since the signals are practically processed in the 
discrete time domain. There exist different ideas on how to increase the 
resolution beyond the sampling time, which can be separated into inter- 
polation techniques in the time domain or the phase domain. While the 
former category is more widespread, it needs a well-chosen interpolation 
model. In contrast, phase-based approaches are built on the assumption 
of narrow-band signals, which means that the measurement signals ap- 
proximately only contain a single frequency. This assumption allows the 
subsample TDE by using linear interpolation or approximation in the 
phase domain to get the phase shifts or the slopes of the phase. 

Let us first consider the interpolation in the time domain. Viola and 
Walker [122] have given an overview of the possible strategies, which in- 
clude interpolation of the signals or interpolation of the cross-correlation 
function. If the signals are interpolated before the cross-correlation, they 
must still be resampled at a higher sampling rate to calculate the cross- 
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correlation. This results in a large computational effort and the estimated 
time delay will still be quantized, though at a finer resolution. Alterna- 
tively, the parameters of the interpolation model can be used to directly 
estimate the time delay, e.g., the root of the analytic derivative of a spline 
interpolation directly leads to an expression for the time delay [122]. In 
the second strategy, the cross-correlation is directly interpolated in the 
immediate neighborhood of the discrete global maximum. The interpo- 
lation models depend on the used pattern-matching function such as 
the cross-correlation, the ASDF, or the zero crossing curve presented 
by Shaswary et al. [112, 113]. In most cases, the interpolations are based 
on parabolic [50], cosine [26], or spline fitting [122], where cosine fitting 
usually performs better than parabolic fitting [26]. Nandi [88] has shown 
that the parabolic fit can also be applied to multiple points of the cross- 
correlation function which offers significant improvements at low SNR. 

For the second group of subsample approximation approaches, the 
phase information is necessary, which can either be extracted by sine 
fitting [95], transformation into the Fourier domain [117] or other model- 
based fittings. For example, Grennberg and Sandell [41] used the cross- 
correlation between the delay signal and the Hilbert transformed ref- 
erence signal to get a direct estimator for the phase. Due to the phase 
ambiguity, only small time delays are allowed if the result is required to 
be unique. Therefore, the phase-based methods are typically combined 
with rough estimators to obtain the integer number of period durations 
contained in the time delay [95, 117]. 

To date, the correlation-based methods and its many variants are con- 
sidered by many researchers as the “gold standard” for TDE [112, 122], 
which is why recent publications still deal with its further improvement, 
such as the frequency-sliding approach presented by Cobos et al. [18]. 


2.1.3 Model-based approaches 


The last group of approaches for the estimation of time delays in acoustic 
measurement systems is based on physical or empirical models for the 
waveforms. Since the model-based approaches are the least common and 
depend heavily on whether a valid signal model can be established, they 
will only be briefly discussed here based on a few selected publications. 
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Yan and Gu [134] proposed a physical model based on solving the 
differential equation of mechanical vibrations. The rising edge envelope 
of the signal model that solves the differential equation is then fitted 
to the measurement signal to estimate the starting time of the recorded 
ultrasonic pulse. A simpler fit of the rising edge is published by McMullen 
etal. [87], where the rising edge was approximated as a parabola, which 
could be fit using two threshold measurements. 

Apart from that, a popular signal model used for the TDE of ultrasonic 
signals is the Gaussian-modulated cosine signal, where the amplification, 
time delay, and width of the Gaussian window as well as the frequency 
and phase of the cosine are the free parameters of the model. Both Gho- 
lami et al. [38] and Lu etal. [75] presented methods to fit multiple over- 
lapping Gaussian modulated cosine signals. While Gholami et al. [38] 
based the fitting on the expectation-maximization algorithm, Lu etal. 
[75] developed a modified Gauss-Newton method. 

In terms of fitting methods, successful applications of the unscented 
Kalman filter [3], the genetic-ant colony optimization [47], or a weighted 
Fourier transform and relaxation-based method [52] have been reported. 


2.2 Transit-time ultrasonic flow meter 


Ultrasonic flow meters, except for some less common types, can be di- 
vided into two main classes: the transit-time UFM and the Doppler UFM. 
Since the latter requires inhomogeneities to reflect the ultrasonic waves 
and is based on estimating the Doppler frequency of the reflected signals, 
it is not further considered in this thesis. In contrast to this, transit-time 
flow meters work for homogeneous single-phase fluids. They are based 
on the effect that an ultrasonic wave sent in upstream direction propa- 
gates slower than in downstream direction. Other UFM variants, such 
as the cross-correlation flow meters, are not considered because their 
special signal processing requirements are only slightly related to the 
topic of this thesis—the accurate TDE for very small time delays in the 
presence of interfering signals. 

The first flow meters based on the transit-time of ultrasonic pulses are 
in use since the early 1950s [118]. Due to the non-existence of moving 
parts, they are characterized by both low energy consumption and low 
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maintenance requirements. Hence, their field of application ranges from 
the food to energy industry, where large throughputs and, therefore, 
large cross-sections are required. While the first UFMs used two paths 
and were built completely analog [34], the systems have been further 
improved to use multiple frequencies [90], phase-shift approaches [137] 
and single path circuits [73]. In order to apply UFM for the measure- 
ment of highly corrosive fluids or in a sanitary environment, so-called 
clamp-on transducers, which allow excitation of the ultrasonic waves 
through the pipe wall without contact to the fluid, were presented by 
Tobin [120] in 1972. Soon after their introduction, Lynnworth reported 
their accuracy limitation due to interfering signals [78]. After years of 
constant improvement, a new chapter was opened in the 1990s with the 
introduction of digital signal processing techniques to the field of TDE 
in ultrasonic applications [32, 92]. Different methods, such as improving 
the robustness of the zero crossing employing the Hilbert transform [99] 
or reconstructing flow profiles using a parametric model [84], were pub- 
lished. With the improved accuracy, the possible applications have also 
further expanded. For example, Schwarz et al. [105] have shown that the 
flow in partially filled drainage pipes can be measured with resolutions 
of less than 1 cms‘. To this end, the UFM system was built to provide 
measurement signals with little noise and interference, and the time- 
delay differences were measured using a high-precision TDE chip that 
provides time-delay resolutions of 55 ps. 

Three different aspects of UFM in which research is still being carried 
out can be named: signal processing, flow profile considerations and 
transducer technologies. 


= In the field of transducers for UFM, Li etal. [72] reported new 
miniature clamp-on transducers, which contain two piezoelectric 
elements, one for the measurement of the time-delay difference 
and one for the measurement of the SoS. 


The influence of flow profiles on the estimation of the average 
VoF was investigated by Wiranata and Ardana [131]. In addition, 
a reconstruction of the flow profile can also be made, as shown 
in the work of Mandard etal. [84]. To this end, a parametric flow 
profile model was fitted to the measurement data from a multipath 
UFM. 
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« The latest research on the signal processing methods featured 
the TDE via improved cross-correlation methods [81], local peak 
fittings [80] similar to the envelope method, and cluster analysis 
of local peaks by processing multiple measurements [126]. 


2.3 Interfering-signals suppression 


As already mentioned in the introduction, the interfering signals limit 
the accuracy of ultrasonic TDEs. If this limited accuracy is not tolerable 
due to the requirements, an interfering-signals suppression must be 
used. In the literature, a distinction is made between algorithmic and 
mechanical suppression where the detailed procedure depends on the 
application. Since the application example in this dissertation is the UFM, 
the SotA suppression approaches for UFM are listed, even though not 
many solutions to this problem have been published. 

In the group of mechanical suppression falls the use of pipe materials 
that are poor at transmitting structure-borne sound, e.g., the plastic pipes 
used in the earliest UFM clamp-on measurements [54]. Apart from that, 
the use of damping mats or other attenuation elements is reported by 
Lynnworth and Liu [79]. 

In contrast, the algorithmic approaches proposed in the literature are 
all based on the fundamental principle that the ultrasonic waves, which 
are transmitted solely through the pipe wall, are stationary, whereas the 
target signals exhibit dynamic behavior because they are dependent on 
the flow [79]. There are several methods to exploit this principle. Roos- 
nek’s approach [99] works by collecting different measurement signals 
whose time delays uniformly cover a period duration to use destructive 
interference for estimation of the stationary interfering signals, i.e., the tar- 
get signals propagating in the fluid cancel each other out when averaged, 
leaving only the interfering signals. However, this approach is highly 
dependent on the availability of the measurement signals, since the time 
delays are required to cover at least one-period duration. Another ap- 
proach has been published by Jacobson et al. [49], where the suppression 
of the non-varying interfering signals is realized by using a high-pass 
with respect to a set of consecutive measurements. Mansfeld et al. [85] 
proposed a similar method, where the signals of consecutive measure- 
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ments are subtracted to remove the stationary signal components—this 
is equivalent to a simple numerical differentiation according to the mea- 
surement index, and thus this method also has a highpass characteristic. 
Nevertheless, the methods proposed by Jacobson etal. [49] and Mans- 
feld etal. [85] require highly fluctuating target signals to get sufficient 
target signal levels, which in turn are necessary to be robust against 
measurement noise. 

In summary, the interfering signals suppression is based on either 
attaching damping materials to the pipes, or using the assumption of 
stationary interfering signals. In the case of the latter, the time delays of 
the target signals must necessarily either follow a uniform distribution 
or change permanently at a fast rate. 
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Since the transit-time UFM is used as an application example in this thesis, 
the fundamentals of the measurement principle, the most important 
process influences and the derived signal models are described in this 
chapter. The focus of this chapter is to present a signal model of UFM 
which describes the classification of the ultrasonic signals into target and 
interfering signals, their sources and properties, and how the influencing 
variables temperature and flow velocity are included. 

Before the measurement system is investigated, an introduction to the 
fundamentals of ultrasonic systems is given in Section 3.1. Based on this, 
the signal model is derived in the following two sections (an overview 
of the sections and their contents is given in Fig. 3.1). In Section 3.2 the 
properties and sources of the target and interfering signals are presented 
and the structure of the signal model is given. Subsequently, the influ- 
ence of the temperature and the VoF on the target and interfering signals 
is investigated in Section 3.3. This finally results in the complete signal 
model with all components, their dependency on the VoF and the temper- 
ature as well as a collection of the properties of the signal components, 
e.g., bandwidth or distribution over the time range. All insights and 
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Figure 3.1 Topics of Chapter 3. 
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assumptions are supported by experimental tests, whose designs and 
implementations are described in Section 3.2.1. 

Finally, the knowledge gained on the properties, influences, and struc- 
ture of the signal model is used in Section3.4 to generate simulated 
data sets. These can be used to validate the intermediate steps of the 
developed filtering methods, which are presented in chapters 4 and 5. 


3.1 Ultrasonic fundamentals 


3.1.1 Ultrasonic flow metering 


The principle of UFM is, as mentioned in the introduction, based on 
the superposition of the SoS c,,; with the VoF v. To this end, two axially 
distanced transducers UT 1 and UT 2 are used to transmit one ultrasonic 
pulse each in the upstream and downstream direction. For this, the 
transducers can be operated both as senders and receivers. Figure 3.2 
depicts a sketch of this principle. The resulting propagation velocity 
is dependent on the path angle a given by the angle between the flow 
vector and the propagation vector. Thereby, the waves in upstream and 
downstream direction propagate with the velocities 


vı = Cy ` v cos(a), (3.1a) 


Up = cy +v: cos(a), (3.1b) 


respectively. Given the path length L, this results in the absolute time 
delays 


L 
T= ————., (3.2a) 
i Cu — v: cos(a) 
L 
Tag = — —— (3.2b) 


Cm +v: cos(a) ` 
By resolving equations (3.2) to the VoF v, we can get the relation 


Taı Ta L 
v= i a, l, (8.3) 


aes 2 cos(a) 
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Figure 3.2 Setup for UFM with necessary system parameters (following [145]). 


which simplifies to 
2 
T ™ 


a NL 3.4 
v È 2 cos(a) a,l Ta,2 ( ) 


if the condition c), > v holds. In liquid fluids such as water with a SoS 
around 1480 ms", this is mostly fulfilled. Otherwise the exact equation 
is to be used, which requires the measurement of the absolute time 
delays and not only the time-delay difference r between upstream and 
downstream direction. 

Either way, only the average velocity of the flow along a single ultra- 
sonic transmission path can be measured this way. To determine the 
mass flow rate, the diameter, the density, and the flow profile have to be 
considered. While the first two are defined by the medium and the pipe, 
the latter depends on the Reynolds number and is usually included in 
the measurement equation via multiplication with a hydraulic correction 
factor. The following flow profile considerations are summarized from 
looss et al. [48]. For fully turbulent flows with Reynolds numbers larger 
than 4000 the flow profile can be described by Prandtl’s law [48] 


v(R) = Umax (1 = aa) (3.5) 


as a function of the maximum velocity in the middle of the pipe v max 
and the distance to the center of the pipe R. Note that this flow profile is 
only valid for a long enough inlet length. The parameter p is dependent 
on the Reynolds number and can be determined experimentally or by 
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an approximated equation. The integration along the diametral acoustic 
path of the ultrasonic wave, whereby diametral means that it passes 
through the center of the pipe, results in an average velocity 

v 


= =., 3.6 
v 1+p (3.6) 


Inserting (3.6) in the integration result over the pipe cross-section 


ge a) (8.7) 


yields the hydraulic correction factor 


E E (3.8) 
v 2 

This correction factor can be used for all diametral acoustic paths. In 

case of multipath UFM, where multiple ultrasonic transducer pairs are 

attached to the pipe, a weighted average of the different single path 

velocities 


_ 1 
v= A 3 WV; (3.9) 


is typically used [84]. Thereby, the weights have to be adapted to the 
path configuration and the expected flow profile [131]. Finally, the pipe 
cross-sectional average velocity can be used to calculate the mass flow 


rate 
D 2 
a=r.:(5) v. (3.10) 


3.1.2 Transducers for ultrasonic waves excitation 


Until now, only the application of ultrasonic waves for flow metering was 
featured. However, the pre-condition of using ultrasonic waves is their 
generation and sensing. To this end, ultrasonic transducers are employed, 
which are mostly based on the piezoelectric effect for high frequencies. 
The exact construction of a transducer depends on the waves to be excited, 
but the core of all transducers is a piezoceramic with electrical connection, 
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backing piezoceramic matching layers piezoceramic damping structure 
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Figure 3.3 Ultrasonic transducer structures for directional longitudinal radiation (left) 
or guided waves excitation, as used for clamp-on UFM (right). The left figure is based on 
[102]. 


backing, and acoustic matching layers. Since a complete listing of all 
transducer types is beyond the scope of this thesis, only one prominent 
type for both longitudinal wave and guided wave excitation is presented 
in this section. 

In Figure 3.3 on the left, a typical structure of an ultrasonic transducer 
to generate pressure waves is depicted. The backing consists of a material 
with a high attenuation that has a similar acoustic impedance as the 
piezoceramic disc to reduce wave reflections. It usually is based on an 
epoxy resin mixed with metal powder [104]. The high attenuation of 
the backing provided by the epoxy resin is necessary to suppress the 
ringing of the transducer which would reduce the bandwidth and affect 
the accuracy, especially for TDE problems. The next few layers are the 
piezoceramic disc and the two electrodes used for electrical connection. 
Finally, several matching layers (for better visualization only two are il- 
lustrated in the figure) with adapted thicknesses are used to transmit the 
ultrasonic waves to the medium with as few reflections as possible [14]. 
However, the above-mentioned relationships of inner reflections, ideal 
layer thicknesses, and attenuations are frequency-dependent functions. 
Therefore, the total behavior of the transducer is highly frequency depen- 
dent. In summary, the size of the layers and their material characteristics 
determine the operational media, the frequency range, and the band- 
width of the transducer. Furthermore, the directivity—consisting of the 
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N 


Figure 3.4 Multipath propagation examples. Depending on the objective to measure, the 
path can be classified into the desired propagation path and the interfering multipath 
propagation. In this scenario, the desired paths are drawn in blue, while the interfering 
paths are drawn in red. 


main and side lobes—is influenced by the diameter and the frequency. In 
this context, a higher frequency or a larger diameter leads to a narrower 
directivity angle [102]. 

The second transducer depicted in Figure 3.3(right) shows a wedge 
transducer that is typically used for clamp-on UFM. Here, the waves 
once again are generated by the mounted piezoceramic and propagate 
as longitudinal waves. After traveling through the wedge material to the 
interface between the transducer and the plate, they get converted to 
shear or Lamb waves propagating in the pipe wall [72, 102]. By designing 
the wedge angles to match the phase velocity of the wedge material only 
to a certain wave mode in the pipe wall, the propagation of other wave 
modes can be reduced. Possibly reflected waves get absorbed by the 
damping structure at the top of the transducer. Subsequently, the shear 
or Lamb waves in the pipe wall radiate into the fluid medium under 
the angle a according to Snell’s law of refraction [23]. In contrast to the 
transducers for longitudinal waves excitation, the wedge transducers are 
mainly used for the generation of guided waves in solid media. 


3.1.3 Multipath propagation 


In many cases, ultrasonic measurement systems are subject to multipath 
propagation, which can be represented by a convolution of the excitation 
signal with a sum of acoustic room impulse responses. This is because 
the excited ultrasonic transducers form side lobes in the directional char- 
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Figure 3.5 System theoretic description of multipath propagation. 


acteristics or the waves can propagate as guided waves in solids, which 
can furthermore lead to a continuous emission of waves. 

Two example scenarios are depicted in Figure 3.4. On the left, a side 
lobe leads to an acoustic path that gets reflected to the receiver by a solid 
object, while on the right, the guided waves propagate in the bottom 
plate, radiate into the liquid medium, and are reflected back. The multi- 
ple acoustic paths obviously have different time delays and waveforms. 
If the objective is liquid level metering, the guided waves propagating 
solely in the plate have to be considered as interfering signals. In Fig- 
ure 3.5, a system theoretic representation of the depicted scenarios is 
shown. Each path is represented by its individual room impulse response 
and the corresponding impulse responses for the sending and receiving 
transducer. 

A more general model of ultrasonic time delay measurements with 
multipath propagations can thus be formulated as follows: 


T 
Salt ) xX (ga alt )* gp alt )) + na(t). (3.11) 


i=0 


In (3.11) the acoustic room impulse response of the system and the sum- 
marized transducer impulse responses for each path are denoted by 
ga i(t) and gr ;(t), respectively, and + denotes the convolution. Note that 
each path has its own transducer impulse response due to the direction- 
specific excitation of ultrasonic waves. If only the time delay of the direct 
path in a homogeneous non-dispersive medium is relevant, the direct 
path impulse response g4 „(t) becomes ö(t — 7) and the total impulse 
response can be separated into 


g(t) = gr(t —7) + Sanat ) * gr a(t), (3.12) 
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where the multipath propagations are modeled by >? gq ;(t) * gp (t). 
Then, the signal components resulting from the convolution of s,,(t) 
with any multipath impulse response can be considered as an interfering 
signal, whose filtering is the main focus of this thesis. 


3.1.4 Lamb waves propagation 


In contrast to fluids, where the predominant mode of propagation is a 
longitudinal wave, solid materials also feature the propagation of shear 
waves [100], where the particle movement is orthogonal to the prop- 
agation direction. Furthermore, the wave propagation is subjected to 
boundary conditions if the solid material is not infinitely extended. A 
particular waveform, that propagates in thin solid plates with thicknesses 
around the magnitude of the wavelength is the Lamb wave, named after 
Horace Lamb, who first presented an analytic description in 1917 [64]. 
Lamb waves can be classified into symmetric and antisymmetric modes 
where in symmetric modes the particles to the left and right of the plate 
center move in phase toward (or away) from the center, and in antisym- 
metric modes, there is a phase shift of 180 degrees. Since the propagation 
of these waves is constrained to solid media, they are a subtype of guided 
waves. Depending on the frequency and thickness of the plate, higher- 
order modes can also occur. 

The analytic solution presented by Lamb can be summarized in the 
characteristic equations that describe the relationship between the angu- 
lar frequency w and the wavenumber k,,. The so-called Rayleigh-Lamb 
equations [100] 


tan(qd/2) 4k "pq . 
e a 72) > TB for symmetric modes and (3.13) 
tan(pd/2) 4k "pq 


i tan(qd/2) 12) -nF for antisymmetric modes (3.14) 
an g qe — Kw 


show that Lamb waves are dispersive, i.e., the phase and group velocities 
are frequency dependent. Here d describes the thickness of the plate, 
while p and q have to be substituted by 


2 2 
p= (=) —k,? and e=(+) =} (3.15) 
EL ET 
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Figure 3.6 Lamb waves group and phase velocities of the first two symmetric and anti- 
symmetric modes in a steel plate with a thickness of 4mm. 


In equations (3.15) the transverse wave speed cy and the longitudinal 
wave speed cz are material constants that are influenced by Young’s 
modulus, the density, and Poisson’s ratio. Thus, the exact dispersion 
characteristic has to be either identified empirically or through precise 
knowledge of the material properties. Note that (3.13) and (3.14) do not 
have an analytical solution and, therefore, have to be solved numerically 
for each frequency w. As the tangent is periodic, the Rayleigh-Lamb 
equations can have multiple solutions for higher frequencies, which 
explains the existence of higher-order modes. Given the solution of the 
relationship between wavenumber and frequency, the phase velocity 


Up(w) = E) (3.16) 
and the group velocity 
ðk, (w) \ t 
Ug) = ( wl ) (3.17) 


can easily be calculated. The resulting phase and group velocities of the 
first two symmetric and antisymmetric modes in a 4 mm stainless steel 
plate with standard material properties are depicted in Figure 3.6. It can 
be seen that the Al and S1 mode can only propagate for f > 0.5 MHz 
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UF 


(a) Test setup for clamp-on UFM systems. (b) Simplified test setup for structural waves. 


Figure 3.7 Setup of the measurement system. 


and f > 0.7 MHz, respectively. Furthermore, the dispersive character of 
the Lamb waves reduces for very high frequencies. 


3.2 Signal model 


The following section introduces the measurement system together with 
its measurement signals and the signal model consisting of the target 
signals and the interfering signals. Measurement signals are recorded in 
different scenarios and examined for their properties such as bandwidth, 
multipath propagations and dispersion. 


3.2.1 Measurement system 


Firstly, the experimental setups used for the identification of the measure- 
ment system are presented. Figure 3.7(a) shows a pair of two transducers 
mounted opposite each other on a steel pipe with a wall thickness of 
4mm and a diameter of 80 mm. Depending on the objective of the in- 
vestigation, water or air is used as the medium. The other setup (see 
Fig. 3.7(b)), used solely for the investigation of structural waves, consists 
of a transducer pair on a 4mm steel plate that is loaded half-sidedly 
with water. Since the plate is from another manufacturer, the material 
is not exactly the same as the pipe wall in the first scenario, but both 
times stainless steel was used. Using the curvature of the plate due to its 
own weight, it can be achieved that only the lower surface of the plate is 
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loaded with water and not the edges or the top surface. This mimics the 
situation in a water-filled pipe without disturbing multipath propaga- 
tions. To this end, the basin is built deep enough and the quadratic plate 
is chosen large enough (~ 1m?) so that all reflections from the bottom 
or the edge of the plate are directly separable in the time domain. The 
transducers are distanced 45 cm apart. In both setups, wedge transducers 
as illustrated in Figure 3.3 were used. 

The measurements were carried out using a PXIe-1062 station with 
a PXIe-5171 ADC module and a PXI-5412 DAC module from National 
Instruments. Both the ADC and DAC modules operate at asampling rate 
of 50 MHz, so all calculations could be performed in the digital domain. 
Unless otherwise specified, the Gaussian modulated cosine signal 


(t —5 ys)” 


t)=10V.- = 
Sy (t) a» ( 2. (1.5ps)2 


) - cos (27 fot) (3.18) 


was used as the electrical excitation, where the center frequency was set 
to 700 kHz. 

Three different data sets were acquired for the system identification 
from the setups presented in Fig. 3.7. For the first data set Dj, 1, the 
setup from Fig.3.7(a) was installed in a water circulation system with 
a pump to control the flow and a reference flow meter to get the VoF, 
which is used to calculate the ground truth for the time-delay difference. 
One measurement each was recorded for the VoFs 0.5ms 1, 0.6m s71 
and 0.7ms!, where each measurement consists of one recording for the 
downstream waves s, (t) and one recording for the upstream waves s ,(t). 
In order to use a uniform terminology for the signals in UFM that is also 
consistent with the SotA methods for TDE, the recordings are called 
delayed and reference signal. As the upstream signal is obviously slower 
than the downstream signal, it is denoted as the delayed signal, whereas 
the downstream signal is denoted as the reference signal. Since the used 
VoFs are quite high and a flow straightener with sufficient inlet length 
was used, a fully developed turbulent flow profile resulted. For the sec- 
ond data set D;, 5, the empty pipe was placed in a temperature chamber 
and while the temperature was increased, 72200 measurement signals 
were continuously recorded, packaged in blocks of 100, and averaged, re- 
sulting in 722 signals. Because there was no water or vacuum in the pipe 
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Figure 3.8 Typical measurement signals of the pipe setup with water (left) and air (right) 
as medium. Only the downstream direction is drawn since both signals are visually the 
same at the present scaling. 


during the experiment, the air is present as the medium. In contrast to the 
first two data sets, the last data set Diq 3 was not recorded from the pipe 
setup, but from the plane plate setup shown in Fig. 3.7(b). Several mea- 
surements were carried out, in which the excitation frequency in (3.18) 
was increased step by step from 300 kHz to 1 MHz with an increment of 
10 kHz. 


3.2.2 Target signals 


All signal components that show the measuring effect—a time-delay dif- 
ference between up- and downstream direction—are considered as target 
signals in this thesis. However, for a better representation of the physical 
interrelationships, the target signals have to be further subdivided based 
on their respective propagation path and transducer impulse response. 
To explain the different physical effects that characterize the UFM system, 
an example signal from the data set Dj, | is depicted on the left in Fig. 3.8. 
In contrast to the scenario where the pipe was filled with air (see right 
diagram in Fig. 3.8), two major signal components starting at about 60 jis 
and 100 ps stand out. Furthermore, other small but widely distributed 
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signal components can also be seen, which already start at 30 ps. A closer 
look at the difference signal 


As(t) = s,(t) — Salt) (3.19) 


between downstream signal s,(t) and upstream signal s(t) at a high 
VoF yields that these small widespread signal components do not contain 
any useful information, since As(t) is nearly zero in their time ranges. 
It is noticeable that the two major wave packets, both of which have a 
significant time delay for high VoF, differ significantly in their envelope 
and especially in their time of arrival. This leads to the conclusion that 
these major wave packets result from a multipath propagation through 
the fluid. As each path has its own path length L; and its own average 
VoF v,, a valid signal model for the target signals can be described by the 
sum 


or ( ma =) l (3.20) 


A CM 


The sign of the time delay 7; depends on whether the signals result from 
the measurement in downstream or upstream direction. It is important 
to note that the exact shape of each signal component s, ;(t) is not con- 
strained through this model, although the individual signal components 
s, i(t) may very well be strongly correlated with each other. The num- 
ber of propagation paths is dependent on the geometric relations of the 
respective pipes and is, therefore, not yet determined here. 

In addition to their propagation behavior, the target signals are charac- 
terized by their bandwidth. That is why the PSD of the recorded signals 
are calculated using Welch’s method [130]. As the window in Welch’s 
method a Hamming window with length 60s for the target signals’ 
PSD and length 10 ps for the noise signals’ PSD was used. To improve 
the noise PSD estimation, the noise signals can be isolated in the time 
domain, since the ultrasonic waves need a minimum travel time, and 
thus the measurement signal contains only noise in the beginning. The 
results are plotted in Fig.3.9. From this diagram, several conclusions can 
be drawn: 
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Figure 3.9 PSDs of the different measurement signals calculated by Welch’s method. The 
PSD of the isolated noise signals (| NV ,.(f) N;(f)|) was calculated using only the time ranges 
that do not contain target or interfering signals. Since the interfering signals are suppressed 
in the difference signal, the PSD of the difference signal (|S A (HS ar f)|) shows the isolated 
PSD of the target signals. Furthermore, the PSD |S,(f) S% (f)| removes the noise influence, 
since n,(t) and nq(t) are uncorrelated. 


= The transducers show a nonlinear behavior to a small extent since 
higher-order harmonics of the excitation frequency are visible in 


ISSN] and |S,(f)Sq(f)] [70]. 


= The noise in the measurement signals can be assumed as white for 
high frequencies because of its flat PSD. The increased values at 
lower frequencies can be explained by either a lowpass character- 
istic of the measurement electronics or residual reverberations of 
ultrasonic waves from previous measurements. 


= Unfortunately, the transducer transfer function limits the possible 
bandwidth, which can be concluded from the fact that the band- 
width of the measurement signals is smaller than the bandwidth 
of the electrical excitation function. 


= Lastly, the noise components in s,(t) and s,(t) are mostly stochas- 
tically independent from each other, which is shown by the dif- 
ference between the PSD |5,(f).S7(f)| and cross power spectral 
density |S,(f)S%(f)|- 
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To summarize, the target signals can be modeled by (3.20), where 
each signal component is individual concerning its waveform s, ;(t), 
its absolute time delay r,; = L,/c,, and its time-delay difference 7;. 
Furthermore, every waveform is limited in its bandwidth. 


3.2.3 Interfering signals 


Apart from the target signals, the measurement signal contains further 
signal components. These remaining signals do not contain the measur- 
ing effect and are considered interfering signals. However, since both 
the target and the interfering signals are always superpositioned in the 
measurement signals, they are difficult to separate using conventional ap- 
proaches. Therefore, a trick is needed to examine the interfering signals 
in isolation. To this end, the fact is used that the attenuation of high- 
frequency ultrasonic waves is significantly greater in air than in liquid 
media [68, 123]. Furthermore, the acoustic coupling reduces through the 
impedance mismatch, further enhancing the suppression of the waves 
propagating through the medium. Using these relationships, the propa- 
gation paths through the medium in the pipe can be almost completely 
suppressed if the pipe is emptied. An example measurement signal re- 
sulting from the data set Diq 9, in which this trick was applied, can be 
seen on the right side of Fig. 3.8. That it has worked can be seen from 
the fact that the two major wave packets seen in D;, , are no longer 
present. Obviously, the present signal is caused solely by waves propa- 
gating through the pipe wall, as this is the only path left when the paths 
through the inside of the pipe are suppressed. Note that the amplitude 
and shape of these structural waves cannot be directly compared to the 
situation with water as the medium because the structural waves are 
not attenuated by the water anymore and the automatic gain control 
amplifies the received waves to take full advantage of the resolution of 
the ADC. 

From the data set Dj, 2 three properties of the interfering signals can 
be deduced: | 


= If the examination of the PSD using Welch’s method is repeated 
for the isolated interfering signals, the same frequency range as 
for the target signals can be observed. Therefore, the presence of 
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interfering signals leads to a spectral interference with the target 
signals. 


= The interfering signal can be considered symmetric, i.e., the signal 
component contained in the upstream signal is identical to the one 
contained in the downstream signal. A quantitative evaluation of 
this property can be obtained by the ratio of the signal energies in 
decibel 


lsat) -s ©) 
10log, JE (3.21) 
: e | 


which yields -40.5 dB for the present data set. 


= There is no time range in which the target signals are the only 
active signals, because the interfering signals are extended over 
the entire time domain. 


The reasons for these properties and a concluding signal model of the 
entire measurement signal, including the target and interfering signals, 
are given in the following four paragraphs. 


Interference with the target signals 


The reason why the interfering signals are in the same frequency range 
is easy to find. Both target and interfering signals are excited by the same 
transducer with the same electrical excitation. This leads to a strong 
correlation between both signal components. Without any nonlinear 
element to provide a frequency shift of the structural waves or the target 
signals, the basic frequency range cannot change. Therefore, the structural 
waves interfere with the target signals. 


Symmetry in up- and downstream direction 


The symmetry of the structural waves follows from the reciprocity of 
the whole transmission system, which can be modeled by a passive, 
linear time-invariant four-terminal network [77]. It is only limited by 
the nonlinear wave propagation in water, the nonlinear behavior of the 
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transducers, and the time-variance of the parameters, since the up- and 
downstream signals cannot be recorded at the same instant. Even though 
these effects are very small, the symmetry is limited to -40.5 dB in the 
recorded data set. 


Time extension of the interfering signals 


The unusually long time extension of the interfering signals is difficult to 
explain, especially since the excitation signal is time-limited to 9 jis. Even 
considering a very long ringing of the transducer, this can only explain a 
time extension in the range of tens of microseconds, but not up to 100 ps. 
To further enhance the understanding of the structural waves, the data 
set Diq 3 was recorded. The dimensions, materials, and setup have been 
chosen to mimic the situation ina water-filled pipe without any multipath 
propagations. Because there is no moving medium and the reciprocity 
property of the acoustic transmission system holds, the direction of the 
wave propagation does not matter. Thus, only one propagation direction 
is recorded and the resulting measurement signals are denoted as s,(t). 
According to Lamb wave theory (see Section 3.1.4), the structural waves 
propagate as Lamb wave modes in the pipe wall, as the wall thickness 
is similar to the wavelength. Therefore, in the time-frequency domain, 
high signal energy is expected at the time-frequency coordinates that fit 
the equation 


t= a , with L=45cm. (3.22) 
valf) 
This was also exploited by Kim et al. [57] for the validation of Lamb wave 
modes. 

The assumption is tested by transferring the measurement signals via 
the short-time Fourier transform (STFT) into the time-frequency domain. 
Although there are other methods to visualize the dispersion character- 
istic, such as a 2D Fourier transform, a dense or at least a sparse spatial 
sampling of the wave field would be needed [46, 138]. In order to cover 
the entire frequency range without amplifying nonlinear effects, multiple 
narrowband excitation signals s,,(t; fọ) with stepwise increasing center 
frequency fo are used instead of only one broadband excitation. In order 
to create an artifact-free detailed time-frequency representation from 
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Figure 3.10 Preliminary study for the qualitative representation of Lamb waves using a 
stacked spectrogram of the plane plate data set Dj. 3. For reference, the group delay at 
the propagation distance 45 cm is drawn for the first two symmetric and antisymmetric 
modes. Left: wide window with high frequency resolution. Right: short window with high 
time resolution. 


the multiple measurement signals, a multi-stage processing is necessary. 
It can be divided into four steps: spectrogram calculation, frequency 
response compensation, normalization, spectrogram fusion. Firstly, the 
absolute square of the STFT, also known as spectrogram, 


2 
S? (t, f; fo) = | fi» *(7 —t)s,(7; fo) e "dr (3.23) 


is calculated, where s,(7; fọ) denotes the recorded measurement sig- 
nal after using an excitation with center frequency fọ. For the STFT a 
Gaussian window w(t) is used to get the best possible time-frequency 
resolution. In the next step, the non-uniform frequency response of the 
transducers is compensated by a weighting function W(f) and subse- 
quently normalized. Fusion of all pre-processed spectrograms is per- 
formed by finding the maximum energy from all spectrograms for each 
discrete time-frequency point. In summary, all obtained spectrograms 
are frequency weighted, normalized, and then fused via the maximum 
operation: 
wi Sst F; fo) 


N max WU) Se Fife) = 
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The resulting time-frequency energy map for two different window 
lengths is depicted in Fig. 3.10. Furthermore, the dispersion relations of 
the first two symmetric and antisymmetric modes translated to (t, f) by 
equation (3.22) are plotted. Since the exact material parameters are not 
known and can, except for the density, not be measured with reasonable 
effort, they are, within the typical limits given in standard tables, chosen 
by an optimization. The energy maxima show good agreement with the 
dispersion relation of all four modes. Now it is also clear that due to the 
highly dispersive S1 and Al mode, the waves can extend very strongly 
in time, even if there is no multipath propagation. Furthermore, if this 
dispersion property is transferred to the situation in a pipe, where a 
multipath propagation is present, the time extension can be even more 
pronounced. 


Final signal model 


In conclusion, only a very general model of the interfering signals s; (t) 
can be assumed, where the symmetric interfering signals are interfering 
with the target signals in the same time and frequency range, but are 
independent from the SoS of the water and the VoF. Thus, the final signal 
model is 


s(t) = >> Sti (t — Tai T;/2) +s) +n, (t), (3.25a) 


salt) = 2, Sti (t Tai T 1/2) + s;(t) + nalt). (3.25b) 


Here, the absolute time delay and the time-delay difference 


L; 
Ti = —, (3.26a) 
iT em 
2v; AT 
To= 5 (3.26b) 
CM 


are determined by the individual path lengths L, and the average VoFs 
along the propagation paths. 
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3.3 Measurement influence of process 
parameters 


In the previous section, the signal model was presented under the as- 
sumption of stationary process conditions. However, ultrasonic waves 
are strongly dependent on the physical properties of the materials, such 
as SOS or attenuation. These properties are in turn dependent on the tem- 
perature, pressure, and other process influencing variables. In particular, 
the VoF and the temperature are the two most important influencing 
variables in UFM. Therefore, to get an estimation about the size of the ex- 
pected effects, the influence of these two parameters on the measurement 
signals is investigated in this section. 


3.3.1 Velocity of flow 


The VoF as the measurand should have a strong impact on the measure- 
ment signals. To this end, the resulting time-delay difference r should 
be very sensitive to variations of the VoF, which, as (3.26b) implies, is 
reached by a large axial distance of the transducers or by a small SoS. 
The axial distance can be influenced, but as the angle at which the waves 
can couple into the fluid from the pipe wall is constrained, the possibility 
of influencing is limited. Furthermore, the SoS depends on the medium, 
which is generally not desired to be affected by the measurement system. 
To sum it up, the sensitivity can only be improved to a limited extent, 
which leads to the point that the impact of the time-delay difference in 
the measurement signals rather than the VoF should be investigated. 
One way to evaluate the time-delay difference without the interfering 
signals is the signal difference since the interfering signals can be consid- 
ered symmetric. Figure 3.11 shows on the left the signal envelope of the 
difference signal, calculated by the Hilbert transform, for three different 
time-delay differences, and on the right the value of three local maxima 
depending on the time-delay difference. In this evaluation, the data set 
Dya,ı was used, where the VoF was stepwise increased with multiple 
recordings per step. Generally, a linear relationship of the envelopes’ 
local maxima is shown, but there is also a considerable variation for each 
VoF step. Another point that stands out is the large differences in the 
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Figure 3.11 Difference between delay and reference signals depending on the time-delay 
difference. 


slopes. This can be explained if the target signals are approximated by a 
narrowband sine signal leading to the difference signal 


As(t) = 8, it + 7;/2 — Tai) — Stilt — 74/2 — Tai) (3.27) 
= Ay + cos(wo(t + 7/2 — T,,;)) 

— A, cos(wo(t — T;/2 — x2) (3.28) 

= 2A, cos (wo(t — to) - sin (wo) 7 (3.29) 


which is valid, if the observed time ranges are small enough (around 
one period). Equation (3.29) proves that the relationship can only be 
considered linear if r, « 1/w, holds. It is also shown that the slope 
depends on both the signal amplitude and the frequency of the ultrasonic 
waves. This also explains, why a higher center frequency leads to an 
improved accuracy of UFM systems. In summary, it can be assumed 
that the information of the time-delay difference is contained in multiple 
propagation paths with some paths being more sensitive than others. 
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3.3.2 Temperature 


The temperature affects the signals in many ways, e.g., the transducer im- 
pulse response, the acoustic coupling or the attenuation in the medium, 
each influence being of different significance. Since the proposed algo- 
rithms in the later chapters do not require complicated temperature 
models to perform reliably, the temperature influence on the signals is 
reduced to a simple temperature-dependent absolute time delay. Simi- 
lar to the target signals, a pulse-specific description of the temperature 
influence for the interfering signals is chosen, leading to the simplified 
signal model 


s,(t;v, T) = > Sti (t = Ta i (T) + r;(v)/2) (3-30a) 
+9 Asiat TaT) nlt), 
I 
sa(t;v,T) = er (-7,,(M —1;(v)/2) (3.30b) 


+ 5 A, Ms, — Ta i(T)) + nalt), 
i 


which incorporates also the VoF, the temperature, and the multipath 
propagations of both the target signals and the interfering signals. In the 
following, the verification of the simplified description of the temperature 
dependence (3.30) is presented. Furthermore, a quantitative evaluation 
of the different time delay variations is given. 


Temperature dependence of target signals 


For the absolute time delay of the target signals, an analytic quadratic 
expression for the SoS in water is reformed using the path length to get 
an estimation for the absolute time delay variation per Kelvin. Using 
one diameter as minimal and three times the diameter as maximal path 
length, the absolute time delay sensitivity 

dr, u L ; deu (T) (3.31) 
dT 


EQ at 


can be enclosed in the interval [55 ns K~*, 347 ns K~'], if T € [19°C, 40 °C]. 
The estimation with the help of an interval was chosen because the ge- 
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ometrical dimensions of the setups vary slightly and only the order of 
magnitude of the sensitivity is relevant. The choice of the temperature 
range is motivated by the process conditions present in the later experi- 
ments. 


Temperature dependence of interfering signals 


In contrast to the target signals, the interfering signals show a much 
smaller temperature dependence, which is more difficult to model. The 
literature reports several methods ranging from simple time-shifting and 
time-scaling [21] to physically modeling the temperature dependence 
of the whole system consisting of the piezoelectric transducer with all 
acoustic couplings and Lamb waves [65]. Although time-scaling has 
proven to be very popular, since it is easy to parameterize and compute 
using the scale transform [43], it cannot handle multipath propagations. 
To overcome this problem, the pulse-specific description of the interfer- 
ing signals in (3.30) is introduced and implemented using the Matching 
Pursuit (MP) algorithm followed by a second constrained MP algorithm 
to track the small change induced by the temperature variation. The 
method was adapted from Lu and Michaels [74], who used it to differen- 
tiate between signal changes due to temperature variation and material 
damages. Since the detailed results of the MP-based modeling of the 
temperature effects have already been published in an own previous 
work [142], the following paragraph only briefly outlines the approach. 

For the investigation of the temperature influence, the measurement 
signals from the data set D;, 5 are used, as they contain isolated inter- 
fering signals subjected to multipath propagation. The first step consists 
of the free decomposition of the measurement signal using the MP al- 
gorithm. The aim is to find a sparse representation of the signal s(t) 
such that the new signal representation compresses the signal energy in 
as few coefficients c, as possible. Coefficients below a specified thresh- 
old are omitted, which equals to them being set to zero. Each signal is 
approximated by a weighted sum of basis functions: 


s(t) = (t) = clk] Y(t). (3.32) 
k 
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The sparse representation problem is formulated by the optimization 


minimize 


K 
s(t) — So clk] v(t) 
k=l 


2 
wrt. d,(t)eG, (3.33) 
2 


where, given a constant K, the best functions Y, (t) have to be chosen 
from an overcomplete function set G = Frame{w, (t), k = —oo,..., o0}. 
The optimization problem (3.33) is also known as sparse approximation, 
which is NP-Hard. Optionally, the problem can be solved by a convex 
relaxation leading to the basis pursuit algorithm or LASSO (least absolute 
shrinkage and selection operator). However, in this case, it is solved by a 
greedy algorithm, iteratively choosing the locally optimal function from 
the frame. The solution may not be globally optimal, but the convergence 
is ensured and the computational effort is reasonable. In summary, the 
MP decomposition of a signal s(t) can be formulated by the iterative 
rules, which are modified from Xu etal. [132]: 


ro(t) = s(t), 


y,(t) = argmax|(r;(t), >, (t))/, 
we lteG (3.34) 


cli] = (rit), Yi), 
Tisi(t) = r,(t) — 2 - Refeli] - ¥,(t)}. 
The decomposition is stopped after K iterations. 
The frame G is built up from energy normalized Gabor wavelets 


(t —t,)? 


2 
20%, 


w(t) = laa” (- 


with variable time delays t,, modulation frequencies f, and time dura- 
tions o}. Note that using the complex Gabor wavelets in combination 
with taking the real part of the projection during the update of the resid- 
ual signal r(t) leads to the optimal phase, which can be calculated by 


Mir) 
vy = arctan (aon a (3.36) 


+ j2n it) 3 (3.35) 
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Figure 3.12 Example residuals of the MP decomposition of interfering signals after differ- 
ent iterations. For better visualization an offset was added to the different residuals. 


The last problem to be solved is how to find the locally best function 
from the frame in each iteration step. This can take up a lot of computa- 
tional effort if a fine discretization of the parameters 0} = |t, fk; 7%] is 
required, since the discretization determines the cardinality of the frame, 
and thus the number of inner products to be computed. A workaround, 
which nevertheless yields a high resolution of the parameter space, is to 
use an optimization method to solve a constrained optimization problem 
instead of a brute force search through all frame functions. This also 
allows for constraints to be set up on the parameter space by adding 
regularization terms. The final optimization problem is defined as 


O; = atg masle) , 6, = [tks Fr; og] , (3.37) 
with the quality function 


J(0,.) = |(ri(t), Y, Ox))| 

gg kr . (fr g fa + ko ` (or u Tan) 

ae ky á (ty = Ven) 

T Rimax 7 (t; _ aa 2 hit), = bax) , 

where h(t) denotes the Heaviside step function, which is used to limit the 


time range in which the decomposition should take place. The strength 
of the regularization can be tuned by the multipliers ky, ko, ky, Ki max 


(3.38) 
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allowing different combinations of constraints. It has to be noted that 
all restrictions can lead to convergence problems. Therefore, the regu- 
larization was used cautiously. An example decomposition of the first 
signal from Dj, 2 can be observed in Fig. 3.12. Here, the process of the 
MP decomposition is shown by plotting the residuals after 0, 7 and 60 
iteration steps. The limitation of the parameter space to the time range 
tmax can be seen very well. In addition, the last residual energy shows, 
that very high compression is possible. 

After the free MP decomposition of one measurement signal, a re- 
stricted MP decomposition of the remaining signals follows. The ob- 
jective of the second constrained decomposition is to track the small 
temperature-induced changes without allowing large jumps in the co- 
efficients. This is important because the MP decomposition does not 
result in orthogonal basis functions leading to large deviations of the 
sparse representation if only the order of the basis functions is changed. 
Therefore, the basis functions and their order resulting from the free MP 
algorithm are fixed and used to repeat the iteration rules (3.34), omitting 
the search for the maximum absolute inner product. Applying this ap- 
proach to the analysis of the signals from the data set D;, o results in 
a set of coefficients c, for each signal, where each signal was recorded 
at a specific temperature. Although the second MP decomposition is re- 
stricted, the residual signal energy is still at —25 dB of the original energy, 
if the temperature range of interest T € [19 °C, 40 °C] is considered (see 
[142] for the model quality at all temperatures). 

The absolute values and the phase values of the largest four coeffi- 
cients c, are plotted against temperature in Fig.3.13. Even though a 
single decomposed signal component cannot fully represent a physical 
propagation path, it is likely that some of the different signal components 
originate from different paths. This is shown by the different slopes of 
the absolute and phase values in Fig.3.13. While the absolute values 
correlate with the coupling, the attenuation, and other effects inside 
the transducer, the phase values indicate the time delay of the compo- 
nent, and thus the SoS. Note that in contrast to the interfering signals, 
no temperature dependence of the absolute values could be observed 
for the propagation paths through the fluid. This leads to two possible 
conclusions: either it is a fictitious effect caused by the superposition of 
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Figure 3.13 Temperature dependency of interfering signals in the MP space for the first 
four basis functions. 


different propagation modes with the subsequent MP decomposition, 
or the attenuation is exclusively due to effects in the Lamb wave propa- 
gation in the pipe wall. Either way, this effect has to be kept in mind for 
the later algorithmic approaches. From Fig. 3.13 the maximum quantita- 
tive effect of the temperature on the absolute value and the time shift is 
shown to be 0.61 % K~! and 11.9 ns K~', respectively. These values are 
calculated by the slopes of the red curve from the left plot and the purple 
curve from the right plot. Since the phase values need to be converted to 
time shifts via their frequency, the excitation frequency of 700 kHz was 
assumed. Comparing the time shift of the interfering signals 11.9 ns K~* 
to the minimum effect on the target signals’ time shift 55 ns K~* shows 
that the assumption of stationarity of the interfering signals can lead to 
estimation errors if the temperature increases. 


3.4 Simulated signals 


In this section, the knowledge gained about the signal model and its 
dependence on process conditions is turned into a simulated data set 
to be used for the development of the filtering algorithms. There are 
advantages and disadvantages of these simulated signals. On one hand, 
the ground truth for the interfering signals and target signals is known. 
On the other hand, the signals are only a simplified version of real mea- 
surement signals and therefore, the performance of an algorithm on the 
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simulated data can only be transferred to real measurement data to a 
limited extent. Nevertheless, using simulated data is a good approach to 
prove the concept of a model-based approach if the assumptions on the 
signal model hold. 

The simulated signals are composed of target signals with three prop- 
agation paths through the fluid, the interfering signals, and artificial 
AWGN. Firstly, a signal form for the target signals has to be designed. 
For this purpose, a Gaussian modulated cosine signal with a center fre- 
quency fọ = 700 kHz and a relative bandwidth of 0.25 is filtered by an 
artificial impulse response of a transducer, which is generated using a 
SPICE-simulated equivalent circuit of a piezoelectric transducer. The 
time delays are set according to Eq. (3.2), with a path angle of 50.12° 
and the three propagation distances L, = 9.36 cm, L = 11.51 cm and 
L, = 14.68 cm. Additionally, different amplitudes are set for the different 
paths to simulate a distance-dependent attenuation. Two types of data 
sets are generated each consisting of 1000 signals: one with a constant 
VoF and one with a rising and a falling ramp for the VoF as seen in 
Fig. 3.14(a). For both types, the temperature is set to start at 19°C and 
to exponentially approach 40 °C. To calculate SoS in (3.2), Lubbers and 
Graaff [76] presented the simple quadratic approximation 
3-1 


cu(T) = 1404.3ms™! + errs 


ms! 
-T — 0.04 rar T: (3.39) 
Although they only claimed a high precision for T € [15 °C, 35°C] a com- 
parison to the data from the NIST webbook [67] shows a good agreement 
in the used temperature range T € [19°C, 40°C]. 

Subsequently to the simulation of the target signals, realistic inter- 
fering signals have to be generated. To this end, a propagating wave 
packet is defined once again by a Gaussian modulated cosine signal with 
fo = 700 kHz. This wave packet at x = 0 is propagated in the Fourier 
domain governed by the dispersion relation 


site) = FHS (fe = 0) -exp(—jk,(w)a)}, (3.40) 


with S,(f,x = 0) denoting the Fourier transform of s;(t) at x = 0. As 
there are different Lamb wave modes, the signal energy of the wave 
packet is distributed to the first two symmetric and antisymmetric modes, 
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Figure 3.14 Simulated data set. 
Table 3.1 Overview of the simulated data sets. 
| VoF constant VoF varying 


s,(t) stationary D D 
s,(t) varying D D 


sim,1 sim,2 


sim,3 sim,4 


then individually propagated by (3.40) and finally superpositioned to 
generate the final realization of the interfering signals. An example signal 
from the generated data sets is drawn in Fig.3.14(b). The three propaga- 
tion paths through the fluid are easily recognizable. Furthermore, the 
interfering signals are extended over the whole time range, as it was also 
observable in real measurement signals. 

In order to test the algorithms presented in this thesis once with 
and without model errors, the interfering signals are simulated once 
as stationary and once with a temperature-dependent amplification of 
0.61 % K~? and time shift of 11.9 ns K~*. These values are chosen to match 
the maximum observed amplification and time shifts in the measurement 
data. This, together with the decision of whether to use a constant VoF 
or a varying VoF, leads to a total of four data sets, which are summarized 
in Table 3.1. 
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on Signal Dynamics 


In this chapter, the first proposed class of algorithms to filter the inter- 
fering signals is presented. They are based on the assumption that the 
interfering signals are stationary, while the target signals show a variable 
time delay due to changes in the SoS induced by the temperature or the 
VoF (see (3.25)). In this thesis, these properties are called signal dynam- 
ics of the target signals. For this reason, the presented methods in this 
chapter will be summarized under the term signal-dynamics method 
(SDM). The first section, describing the requirements and the prelimi- 
naries, contains a new representation of the signal dynamics in form of 
point clouds. Subsequently, two methods to process the resulting point 
clouds are introduced in the following two sections leading to several 
methods for the estimation of the interfering signals. Both methods only 
evaluate the direct propagation path. Finally, each method is tested on 
the simulated data sets and concluded by a discussion of the advantages 
and disadvantages of the respective method. 


4.1 Point cloud representation of consecutive 
measurements 


4.1.1 Requirements and preliminaries 


The precondition to use signal dynamics for filtering stationary inter- 
fering signals is to have a package of multiple measurement signals. 
Furthermore, this package needs to be recorded while the varying pro- 
cess conditions lead to time-shifted target signals. In summary, let there 
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be M measurements each for the reference s, (t) and delayed signal s q(t) 
combined to the vectors 


s,(t) = Is. (t,1), ...,s.(t, m), ... ,s,(t, M)]F, (4.1a) 
salt) = [sq(t, 1),...,8g(t,m),..., Sa(t, M)]', (4.1b) 


where the measurement index m denotes a single signal within the pack- 
age. In the following, the process conditions T, v or the derived variables 
T,,7 during a single measurement of the package are always denoted 
by the measurement index, e.g., Tm, Vm» Ta,m: Since the processing takes 
place in the digital domain, the signals are only available in the time- 
discrete form 


S, = (s,(tn,m)) 
Sa = (salta m)) ERM, (4.2b) 


mn 


ge RAN, (4.2a) 


with t,, = nt, being the discrete time sampled at the sampling rate 1/t,. 
To highlight the difference between quantities with physical dimensions 
such as t,, and integer quantities such as n, round brackets are used in 
the physical case and square brackets in the integer case. Furthermore, 
multiple equal independent variables that influence a signal or function 
are separated by comma in the argument, whereas hyperparameters or 
process conditions are separated by semicolon. 

The signal model used for SDMs is simplified to contain only a single 
propagation path of the target signals, i.e., the sum in (3.25) is omitted. 
Cutting out only a small time range, in which the direct propagation 
path is included, is necessary to improve the validity of this model. To 
this end, a preprocessing searches local maxima in the envelope of the 
measurement signals and their difference signals As (t). By introducing 
several regularizations in combination with choosing the local maximum 
with the smallest transit time leads to a robust detection of the correct 
time range. The regularizations used are peak prominence, peak height, 
necking in the envelope, and physically possible limits. The remainder 
of this chapter assumes that the time range, which the SDMs have to 
process, is known or determined by a preceding method. 
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Figure 4.1 Hilbert transform example and 2D representation of an analytic signal. 


4.1.2 Fundamentals on the Hilbert transform 


Real-valued signals are characterized by the fact that the absolute values 
of their Fourier transform are symmetrical around f = 0. In contrast, 
the Fourier transform of analytic signals is non-zero only for f > 0. 
For sinusoidal signals, this has multiple advantages, such as an easy 
determination of the phase, the instantaneous frequency, or the instanta- 
neous amplitude. For these reasons, the measurement signals are firstly 
transformed into their analytic counterparts using the Hilbert transform 
[53] 


Sas(t) = s(t) +j-H{s(t)} (4.3) 
ee, wie Aa: (4.4) 


nt 


However, since the signals are time-discrete, a discrete calculation of 
the Hilbert transform is necessary. This can be realized by calculating 
the discrete Fourier transform (DFT) of the signals, replacing the DFT 
coefficients that correspond to negative frequencies with zeros and finally 
transforming the resulting DFT values back to the discrete time domain 
by using the inverse DFT [86]. 

The resulting analytic signal is complex-valued as (4.3) implies. There- 
fore, in order to visualize the signal, it can be displayed either separately 
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by real and imaginary part or in a 2D area where the real and imaginary 
part are the axes. Figure 4.1 shows an example of the Hilbert transform 
applied to a signal from the simulated data set Dim 1- Itis easy to see that 
the Hilbert transform produces a 90-degree shift of s(t,,), which leads 
to spiral trajectories in the complex area. Since the sampled values at 
each point in time are projected onto their real and imaginary values, the 
direct relation to the time vanishes in the complex area, which makes this 
representation invariant to time scalings or time shifts. The formulation 


Sag(t) = Alt) Pr) (4.5) 


shows that the phase position and the envelope can be calculated taking 
the argument and the absolute value of the analytic signal. Note that if 
the envelope is constant, the trajectory will form an ideal circle. 


4.1.3 Point clouds 


Applying the time-discrete Hilbert transform, introduced in the previ- 
ous section, to each measurement signal from the given signal package 
(4.2) leads to the complex-valued matrices S, as; Saas € CMXN Subse- 
quently, to interpret the signal package as point clouds, the real and 
imaginary parts are considered as abscissa (x-coordinate) and ordinate 
(y-coordinate), respectively, which yields one 2D point cloud per discrete 
time index n for both the reference and delay signals 


Re els! . 
e oag | ene, (4.6a) 
Refe!sT 

e{en d,as | € RM, (4.6b) 


P.[n] = | 


Pain] = Im{eTsT 


a ~ d,as 


where e,, denotes the n-th unit vector and each point cloud consists of 
M points. 

An example of the point cloud representation is shown in Fig. 4.2. If 
the signals and their Hilbert transforms on the left side are sampled 
at the four discrete time steps indicated by the dashed lines, the four 
point clouds depicted on the right side are obtained. Each point cloud 
consists of as many points as there are measurement signals. The length, 
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Figure 4.2 Point cloud generation from measurement signals package. 


curvature, and orientation of the point clouds is determined by the signal 
dynamics, the frequency, the amplitude, and the instantaneous phase 
at the corresponding time step. Although a high SNR was set in the 
simulated data sets, the influence of measurement noise is shown by the 
slight deviation of the point clouds from a perfect line. 

Lastly, the origin of the effect that the point clouds of consecutive time 
steps t,, are misaligned is discussed. The following considerations apply 
to both the reference and the delay signals. If the signal model (3.25) with 
stationary interfering signals and varying time delays 7,, is assumed, 
the signal package (4.1) can be reformulated to 


8,(tn) = [si(tn eh see St (tn 7" + Ian ` si(tn) +n, (tn) (4.7) 


In (4.7), the M-by-1 column vector filled with ones 1,7 = {1} indicates 
that the interfering signals s;(t,,) are identical in each measurement of 
the package. Since the effect is the same, it does not matter if the time 
delay 7,,, is caused by a varying absolute time delay T, m, a varying time- 
delay difference r,, or a combination of both. Repeating the necessary 
steps to get point clouds for this model results in 


P,[n] = P,[n] + piln]: In + Pa, [n] (4.8) 


This shows that each point cloud P,[n] is composed of a shape influ- 
encing component P,;[n], an offset p;[n] determined by the interfering 
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signals, and a noise-induced stochastic component P „ [n]. For further 
investigation of the point cloud shape, suppose that the time delays 7,,, 
in (4.7) cover the interval [t),,¢,,,]. Thus, the component P,[n] can be 
considered as a discrete sampling of the analytic target signal in the 
interval t,, = [t,, — tub ; tn — typ]. For this reason, the shape of P,[n] 
is determined by the shape of the analytic target signal in the interval 
t „. Transferring this concept to other time steps, the next point cloud 
P,[n + 1] shows the shape of the analytic target signal in the interval 
tia = [nai — tub: tri — typ] and so on. A sufficiently small sampling 
time or a sufficiently large interval 7 € [ty,,t,,] leads to a non-empty 
intersection t = t,, N t,,,,. In other words, the shapes of consecutive 
point clouds overlap in cases where the offset p;[n] is zero. Otherwise, 
the misalignment 


piin] = piin + 1] — pin] (4.9) 
is present between the point clouds P,[n] and P,[n + 1]. This concludes 
all effects that are observable from the example point clouds. 


4.2 PCA-based point cloud processing 


This section describes how the point clouds presented in the previous sec- 
tion are processed by the PCA to filter the interfering signals in TDE. The 
point clouds can be processed in two ways: either they are used to first 
geometrically estimate the interfering signals followed by a simple sub- 
traction to improve conventional TDE methods, or the TDE is estimated 
directly from the point clouds. Both methods assume a narrowband sig- 
nal model with limited bandwidth, which was shown in Section 3.2.2 to 
be in good agreement with the measurement signal properties. 


4.2.1 Geometric estimation of interfering signals 


The theoretical motivation to estimate the interfering signals geometri- 
cally is based on the narrowband signal model in analytic signals form 


Sa T) = Ay TH +5, W(t) Anal), (4.10a) 
Bl T) = A, erfor?) +8; as(t) T Na as(t) , (4.10b) 
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with the stationary interfering signals s; a(t). In (4.10), the absolute time 
delay was omitted and a constant frequency was assumed, but neverthe- 
less, the following argumentation for the interfering signals estimation 
is still valid, since the only precondition is the presence of the signals’ 
derivative at two distinct time delays. In order to make the calculations 
more manageable, the analytic signal components are transferred to the 
2D area equivalent to the point cloud representation. To this end the 


operation 
rGHR gG ee (4.11) 


is introduced, which maps a complex number to a 2-by-1 vector. Applying 
the vec operation to all components from (4.10) 


P,(t;7) = vec (s,ac(t57)) ,  Palts7) = vec (Sa.as(t57)) » 
pı(t) = vec (A, e7), p(t) = vec (si as(t)) 


Pn, (t) = vec (ne aslt)) 3 Pr, (t) = vec (naaslt)) (4.12) 
results in the 2D signal model representation 

P(t; T) = pelt +7/2) + Pill) + Pa, (t), (4.13a) 

Palt; T) = Pe(t— 7/2) + pit) + Pn, (t) - (4.13b) 


Remembering the point cloud representation of the signal package, 
short sections of the target signals’ trajectory are given by the point clouds, 
where the lengths of the sections are determined by the signal dynamics 
contained in the package. Even if the true trajectory of the target signals 
cannot be restored due to the sections being too short, it is possible to 
calculate the tangents to the observable trajectory pieces. These tangents 
v,, vg are proportional to the derivatives of the signals with respect 
to the time delay 7. By inserting the complex exponential function in 
the derivative and using Euler’s formula, a relationship between the 
derivatives and the target signals can be established: 


d i = 

Tedin = Tho ; | Pilt + 7/2) x v,, (4.14a) 
dp,(t;r 

un =Tfo ee A -p,(t — 7/2) x va. (4.14b) 
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The rotation matrices in (4.14) are caused by the multiplication of the 
target signals with the imaginary units, resulting from the differentia- 
tion. Subsequently, by inverting the rotation matrices, the normals to the 
tangents 


pilt + 7/2) =c b i “Vv, =: cv, (4.15a) 
0 —1 \ 
pilt — 7/2) = ca ki 0 | "Va =! —Ca V] (4.15b) 


are given for both distinctive time delays 7/2 and —r/2. Note that for a 
better overview, the time dependencies of the tangents v,, Vq and the 
normals vý, v are not displayed. Finally, an estimation of the interfer- 
ing signals can be obtained by inserting the normals (4.15) into (4.13), 
neglecting the measurement noise, and reformulating the equation to 
p;(t). The resulting estimation is given by 


Bl) = p,(t;7) — [1,0] -[ve, va]: (P(t; T) — Palt; 7)) vr, (4.16) 


which leads to the estimated interfering signals §,(t) = e| p(t). Taking 
into account the measurement noise, the estimation deviates from the 
true interfering signals values according to 


B,(t) = pilt) + Ph, (t) 
+ [1,0]-[vr, val? - (pn, — Pn, (4) vi. (4.17) 


From this equation, two conclusions can be drawn: 


= Due to the multiplication with the normals matrix and the projec- 
tion onto the normal v, the x- and y-coordinate of the estimation 
error are not uncorrelated anymore. 


= The variance of the noise gets amplified if the normals matrix 
determinant tends to zero, i.e., the normals v} and vj need to 
have distinct angles, which can be improved if the distinct time 
delays 7/2, —r/2 are further apart. 


This concludes the time-continuous formulation of an estimator for the 
interfering signals based on the tangents to the trajectories. 
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The following is the practical implementation of the estimation based 
on the given discrete-time signals and their derived point clouds P [n], 
P a[n]. Instead of solving (4.16) for all times t, the interfering signals need 
to be only estimated for the discrete time steps n. To further reduce the 
effort, the time range can be limited to only the time range of the target 
signals, where the SNR is higher than a predefined threshold. After 
determining for which time steps (4.16) is to be solved, the next task is to 
calculate the tangents. For this purpose, the point clouds are processed 
using the PCA. The PCA adapted to the processing of 2D point clouds 
can be summarized in two steps. At first, the point clouds at the discrete 
time steps are processed via the rule 


P[n] = P[n]- (tu = Fav) , (4.18) 


with I,, denoting the M-by-M identity matrix, for the reference and 

delay signals respectively, which leads to the offset-free point clouds 

P, [n] and P a[n]. In the second step, the eigenvalue decompositions 
Din] = V,[n]-A,[n]- Vy [nr], (4.19a) 
Dafn] = Valn]- Auln] - V2 [n) (4.19b) 


of the estimated covariance matrices 


5, [n] = = P [n]: P’fn] e R??, (4.20a) 
3, [n] = 5 B,In]- Pin] e R?? (4.20b) 


are calculated. For further information on the PCA and its multiple 
applications see Beyerer etal. [7]. In the present case, the eigenvector 
matrices V,[n], Va[n] are composed of two eigenvectors v € R? each. 
The respective eigenvector with the largest corresponding eigenvalue, 
i.e., the largest corresponding diagonal entry of A, is called the first 
principal component and signifies the direction in the 2D area with the 
largest variance of the point cloud. Therefore, the resulting matrices 
V,.[n] and V ,[n] can be used to get an estimate for the tangents by using 
the respective first principal components v,[n] and vg[n]. To complete the 
necessary variables for the time-discrete variant of (4.16), only the vectors 


59 


4 Novel Separation Methods Based on Signal Dynamics 


, Um{Sas(tn)} xo Ser silta) k 
Be 5,(tn) 
k 0.05 = ge 
0.4 + > \ 
= 0.00 | 
= 
0.2 + P J 
—0.05 |- + 
Re { sS as (t n ) } | | 
H - > 
0.2 0.4 70 72 74 
t„/ns 
(a) Geometric estimation of interfering sig- (b) Interfering signals estimation using the 
nals at t,, using a constant amplitude as- PCA-based geometric estimation with 
sumption. (5,) and without (8f°™) constant ampli- 


tude assumption for the data set D 
Ground truth is drawn in dashed. 


sim,1' 


Figure 4.3 Example of the PCA-based interfering signals estimation. 


p,[n;7] and paln; 7] are missing. These can be approximated using the 
mean value of the respective point clouds 


1 

p,[n;7] = welt ‘Inn, (4.21a) 
1 

Paln; 7] = Mr al] Inn: (4.21b) 


An example of the approach for a single time step is shown in Fig. 4.3(a). 
There are four components shown: the interfering signals p;[n] (green 
arrow), the target signals of the reference and delay signals p,, pg (blue 
and red arrows), their respective point cloud (blue and red points), and 
the resulting first principal components (dashed lines). This figure also 
explains why the solution of (4.16) can be interpreted as a geometric 
intersection between two normal vectors. Repeating the estimation for 
all time steps t,, € [70 us, 75 us] leads to the estimated interfering signals 
5,(t„) in Fig. 4.3(b). However, it is easily observable that the deviation 
from the ground truth is not to be neglected. There are several sources 
for the estimation errors: the measurement noise, non-uniform sampled 


60 


4.2 PCA-based point cloud processing 


point clouds, the tangent estimation, model errors such as the violation 
of the constant amplitude assumption. 

Since the assumption of constant amplitude is the most restrictive 
condition, it is relaxed by adjusting the normals estimation. The proposed 
normals correction is presented in the following paragraph, where due 
to the redundancy of the equation only the equations for the reference 
signals are displayed. The delay signals are processed equivalently. Again, 
the derivation of the adjusted method is based on the signal model. 
Adding a non-constant amplitude to the signal model (4.10) leads to the 
extended model 

Sy aslt; T) = A(t Gi T/2) er folt+r/2) +5; as(t) a. Ny as(t) : (4.22a) 
It follows that the derivative of the signals 
dsr aslti T dA(t + 7/2 ae 
nn = +jA(t + r/2)r fo) ei? Folt+7/2) (4.23) 
T F 


is extended by the derivative of the amplitude, whereby the rotation by 


90° in (4.14) becomes a rotation by the yet unknown angle y leading to 
dp.(t;r) _ [eos(p) —sin(y) 
Zr hae cone pilt + 7/2) x v,. (4.24) 


Therefore, the definition of the normal vectors is changed to become the 
multiplication of the tangent p, with the inverse rotation matrix, i.e. 


ee Be = vy, = evi. (4.25) 
From (4.23) the angle 
u 2nfoAlt + T/2) 
p=arctan a) (4.26) 


can be deduced. However, since only the whole measurement signal with 
all additive components can be observed, the amplitude is not known 
before the interfering signals are subtracted. This can be solved by using 
the signal difference between the analytic reference and delay signals. 
The absolute value of this so calculated difference signal 


[Asas (ti P| = I82,as(&; T) — Saas(ts7)| © 2| sin(afo7)| AC) 


(4.27) 


r,as 
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is an approximation of the time-dependent amplitude curve, if the time 
delay 7 is very small compared to the time duration of the envelope 
A(t). Even though the absolute value of the amplitude curve remains 
unknown, for the estimation of the angle y only the ratio is needed. 
Therefore, the angle is estimated by 


27 fo|As asht; T)| ) 


IAs P/d eee 


9 = arctan ( 


where the derivation is performed in the discrete time by a first-order 
divided difference. 

In summary, for a non-constant amplitude, the approach is nearly the 
same as using the constant amplitude assumption. The only difference 
lies in the calculation of the normals, where an adapted angle ¢ is in- 
troduced. The result of this relaxation to non-constant amplitudes is 
displayed in Fig. 4.3(b) as s*°""(t,, ). The expected improvement of the 
adjusted estimation is proven, as the deviation from the ground truth is 
significantly smaller than when assuming a constant amplitude. 


4.2.2 Interference-invariant time delay estimation 


Apart from the method of first estimating the interfering signals and then 
simply subtracting them, the time delay between reference and delay 
signals can also be calculated directly from the point clouds. To this end, 
the stationary interfering signals are separated from the target signal by 
transforming the measurement signals into a space that is invariant to 
stationary signal components. Since this approach is already published 
in [143], a brief description of this transformation and the subsequent 
post-processing to determine the time-delay difference 7 follows. 
Stationary components in the point clouds have only an influence 
on the offset as proven by (4.8) in Section 4.1.3 and therefore have no 
influence on the shape or orientation of the point clouds. It follows that 
the offset-free point clouds retain only the shape influencing components 
and the orientation. Thus, by using the orientation of the offset-free point 
clouds, expressed by an angle, the interfering signals can be suppressed. 
However, since the measurement noise is also still included, a reliable 
estimate can only be obtained if the time delay variation of either the 
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Figure 4.4 PCA-based point cloud processing to get interference-invariant phase values 
of an estimation of the interfering signal. 


absolute time delay r, or the time-delay difference rin the signal package 
is significantly larger than the variance of the measurement noise. 

In order to calculate the orientation of the point clouds for both the 
reference and delay signals, the PCA is employed. From the respective 
first principal components, given by the eigenvectors v,[n], vain] € R? 
with the respective largest associated eigenvalues, the orientations are 
determined by 


T . 
pcan] = arctan (See) j (4.29a) 
r [n]: 1: 
T . 
®pca,aln] = arctan (oe) (4.29b) 
valn] -e 


After solving the PCA equations (4.18)-(4.20) and calculating the an- 
gle (4.29) for each time step n, the measurement signal is completely 
transformed into a space invariant with respect to stationary signals. 

However, the direction of each eigenvector is ambiguous, i.e., it is 
not possible to influence the PCA whether the eigenvector is rotated by 
180 degrees or not. Due to this problem, there are still random phase 
jumps present in the phases calculated according to (4.29), which need 
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to be removed before the phase values can be further processed. For this 
reason, a case distinction is included in the phase calculation, adding 
r if v'[n] - e, < 0. Subsequently, the phases thus obtained are mapped 
linearly to the range [—7, 7]. In summary, for both the reference and delay 
signals, the phase calculations in (4.29) are adapted to 


in =r v,[n] > 0, vo[n] > 0 
épca [n] = 2- arctan Gar +4 m v,[n] >0, vn] <0 , (4.30) 
0 v,[n] <0 


where the components of the eigenvectors v"[n] - e} and v"[n] -e, are 
abbreviated to v,[n] and v,|n], respectively. 

The course of the orientations is shown in Fig. 4.4(a). It should be noted 
that the orientations contain only a reliable information in time ranges 
with a significant level of target signals. Outside these time ranges, these 
orientations are determined only by the measurement noise and carry 
no information about the time-delay difference. The time range shown 
in Fig. 4.4(a) was selected on the basis of the transit time of the direct 
propagation path. Furthermore, the phase values obtained by directly 
using the Hilbert transform are included in the plot of the phase values. 
It is observable that the solution for the ambiguity of the eigenvector 
direction leads to double slope compared to the Hilbert phase ¢,, [n]. 

To extract the time-delay difference r from the ¢pc, In], dpca all 
orientations, they are separated by the procedure proposed by Kupnik 
etal. [63] at the 27 jump points and approximated by linear regressions 
Ppoa rln t; Open aln t,), as shown in Fig. 4.4(b). The best linear regres- 
sion is selected based on the eigenvalue ratio of the associated orienta- 
tions. This ensures that the linear regression with the best SNR is used. 
The intersection points with the time axis can then be determined from 
the linear regressions, and thus, by simply taking the difference, the 
time-delay difference is calculated. 

However, a closer look at the method shows that only one time-delay 
difference per point cloud can be estimated, i.e., only the average time- 
delay difference of a whole signal package results from this method. 
Further information and a quantitative evaluation of the performance 
can be found in [143]. 
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4.2.3 Discussion 


Both PCA-based processing methods are tested on the simulated data 
sets. Before they can be applied, compliance with the preconditions of the 
SDMs has to be ensured, i.e., the measurement signals need to be sepa- 
rated into signal packages. The included signal dynamics per package is a 
free to choose hyperparameter that controls the robustness against noise, 
but also limits the precision of the tangent estimation since the PCA is 
designed for linearly shaped point clouds. Furthermore, the dynamic 
behavior of a measurement system based on this concept is limited, since 
the interfering signals can only be reliably filtered after the required num- 
ber of measurements has been recorded. For the proof of concept, the 
signal dynamics per package is set to 50 ns, which corresponds to a phase 
shift of 12.6° at a center frequency of about 700 kHz. To obtain estimates 
of the time delays for all measurements, the signal packages are selected 
with overlap such that the entire range of measurements was uniformly 
covered. Because the geometric approach estimates and subtracts the 
interfering signals, before the zero-crossing-based TDE approach from 
Kupnik et al. [63] is applied, this SDM enables a TDE for each signal 
within the package. On the other hand, the interference-invariant SDM 
yields only one TDE per package. Due to the overlapping signal packages, 
the geometric SDM results in multiple TDEs for the same measurement, 
which are therefore averaged to get comparable TDE curves. For easier 
notation later in Chapter 6, the set of hyperparameters for the PCA-based 
SDM 


gPCA = {Ar, ConstAmp} (4.31) 


is defined as the collection of a signal dynamics per package parameter 
Ar and, in case of the geometric approach, a Boolean value ConstAmp, 
which indicates whether or not the constant-amplitude assumption is 
used. 

The results of the proposed PCA-based SDMs are pictured in Fig. 4.5. 
As already expected, the interference-invariant approach performs badly 
for data sets with fastly varying VoFs, e.g., D Dim, Since it can only 
calculate one time delay per package, which limits the allowed dynamic 
behavior. However, if the stationary interfering signals assumption holds 
and the rate of change is slow enough, this approach is most precise and 


sim,2/ 
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Figure 4.5 Results of TDEs using the PCA-based SDM for the simulated data sets 
D D 
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most robust against noise. The geometric approach can alleviate the lim- 
ited dynamics, which is proven by the fact that the results also perform 
well in data sets with varying VoF. However, the expected improvement 
using the non-constant amplitude extension did not occur. One possible 
explanation is that the incorrectly estimated interfering signals when 
assuming a constant amplitude are phased in such a way that they have 
only a small negative effect on the TDE. It is actually also the case that 
the extended method often shows worse behavior both against noise and 
against the systematic error. While the behavior against noise is to be 
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expected, since the estimated correction angle y depends on the noisy 
difference signal, the systematic error can only be explained by the phas- 
ing of the residual interfering signals after subtraction. A non-oscillating 
constant positive TDE error indicates that the residual interfering signals 
are after correction always in 180 degrees phase shift to the target signals. 
This leads to the explanation that the amplitude estimation used for the 
angle determination is the source of error, since the difference signal is an 
input quantity moving with the process influence quantities to the same 
extent as the target signals. The last insight that can be gained from the 
figure is the quantitative influence of non-stationary interfering signals. 
If the interfering-signals variations—with the quantitative descriptions 
determined in Chapter 3—are assumed, the TDE error is within ~ +2 % 
(see lower left of Fig. 4.5). 


4.3 B-spline-based point cloud processing 


The performance of PCA-based point cloud processing depends on the 
extent to which the point clouds are linearly shaped. Therefore, the 
maximal signal dynamics per package have to be limited, which can be 
achieved by forming packages in which the phase of the target signals 
changes only by a maximum of 25°. However, using higher signal dy- 
namics would allow better robustness against noise. Since the shape of 
the point clouds becomes more similar to sections of spirals, if the signal 
dynamics are increased, a curved approximation becomes necessary. In 
this thesis, this problem is solved by applying B-splines to the approxi- 
mation problem. They are fast to calculate and can approximate every 
multidimensional line as a parametric curve where the continuity of the 
curve and its derivatives can be easily determined by manipulation of 
the knot vector. Furthermore, the application of B-splines removes the 
requirement of having uniformly distributed points within the point 
clouds, which in fact only depends on the process conditions during the 
available measurement signals. In this section, the B-spline-based SDM 
that was briefly published in [141] is described in detail. The presenta- 
tion of the approach is structured as follows. Firstly, the fundamentals of 
B-splines are introduced, followed by two 2D approximation approaches 
to get robustly estimated B-splines. Subsequently, one algorithm to es- 
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timate the interfering signals and one algorithm to estimate the time 
delays directly as in Section 4.2.2 are presented. Finally, the concept is 
proven and discussed by application to the data sets D Dim, Where 
the signal model holds, and also tested against signal model errors by 
application to D,,,,3 and D 


sim,17 
sim,4' 


4.3.1 B-spline fundamentals 


The following fundamentals on B-splines are a summary of selected 
contents from [25]. B-splines, also known as basis splines, are piecewise 
defined polynomial functions with limited support, whose derivatives 
have adjustable continuity properties. Here, basis functions mean that 
any other spline function can be expressed as a linear combination of 
B-splines for the same polynomial order. A unique set of B-splines is 
given by the B-spline order l and a knot vector ry = [ro,..-,"7_ı]" with 
Yj41 2 r; defining the limits of the polynomial pieces. Using the given 
order and knot vector, the B-splines can be defined by the Cox-de Boor 
recursion formula [20, 25] 


1 or 
balr) = re Pate) , (4.32a) 
0 else 
rer, Tiy TT 
bi pal) = —— h, (r) ns (4.32b) 
one Tea TG Teste Fir seas 


The resulting set of functions b; ;(r) are called B-splines and fulfill the 
sum-to-one property 


Solr) =1, Vr E€ [ro rr) (4.33) 


within the interval [r,,r 7—1) spanned by the knot vector r,,. Furthermore, 
from the recursion formula (4.32) it can be seen that the number of 
intervals [r;,7;,,) in which a B-spline is non-zero depends on the order. 
At order l = 1 this is only one interval, and then the number of intervals 
increases with the B-spline order ! at the same rate. This means that, 
for example, B-splines consisting of quadratic polynomials (l = 3) are 
spread over three intervals. From the same relation, it follows directly 
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that for a B-spline of order l at least l + 1 knots have to be present, which 
enclose / intervals. To ensure that this is always the case, the first and the 
last element are usually appended l — 1 times to the given knot vector 
resulting in 


1-1 I 1-1 
— — [L m ——— 
Py = Ir. ToTo Ties Tia iate (4.34) 
mm — U 


i+1 141 


It should be noted that successive knots may be identical (r; = r;}1), Le., 
the intervals can be empty. In summary, the properties of the B-splines 
and its derivatives can be derived as follows based on the given knot 
vector: 


= In non-empty intervals, the B-splines are polynomials of order 
i—1. 


If the end knots ro, r;_, of the given knot vector are unique and 
appended l — 1 times at both sides as in (4.34), there are a total 
of F = I + l — 2 B-splines, resulting in F degrees of freedom for 
subsequent interpolation or approximation tasks. 


= The B-splines are continuously differentiable at the knots up to the 
(l — 2)-th derivative, if the knots defining the B-spline are distinct. 
For each repeated knot, the number of possible derivatives at the 
corresponding value decreases by one. 


An example of quadratic B-splines (l = 3) is shown in Fig. 4.6(a). It is 
observable that in each interval three B-splines are non-zero. Due to the 
given number of knots and the B-spline order, there are six different 
B-splines present in this example (F = T + 1 — 2). Except for the end 
points, the first derivatives of the B-splines are continuous at the knots, 
i.e., atr € {1,2,3}. Additionally, this example illustrates the necessity to 
add knots at the beginning and end of the knot vector, because without 
these additional knots, the splines in the intervals [0, 1] and [3, 4] do not 
have enough degrees of freedom to meet possible boundary conditions 
such as slope or endpoints. 
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(a) Quadratic B-splines with knot vector (b) 2D approximation of points. 
ry = (0,0, 0,1, 2,3, 4, 4, 4]. 


Figure 4.6 B-splines and application to an approximation of 2D points. 


Finally, the B-splines can be used to generate any 2D spline curve 
through a linear combination 


F-1 


ra 2 Hi bil), (4.35) 


i=0 Lyi 


with the control points [c,, ;, Gral Of course, curves can also be gener- 
ated in only one or more dimensions, but since in this work B-splines 
are used to approximate 2D point clouds, only the representation of 2D 
spline curves will be discussed here. The results of applying B-splines to 
approximate 2D points are depicted in Fig.4.6(b). Note that this curve 
is only possible due to the representation as parametric curve, because 
a conventional function y(x) can only map a single y-value to a given 
x-value. 

Another advantage of B-splines lies in the possibility to calculate the 
derivatives analytically. Since the B-splines are defined recursively by 
(4.32), it follows that the derivatives of B-splines can be calculated using 
the B-splines of order ! — 1: 


db; (r) bi (r) bi+1,1-1(R) 
u ee ec 


(4.36) 


Tiy Ti Viste Fir 
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4.3.2 2D B-spline approximation 


Using the B-spline theory, the shapes of the point clouds calculated from 
the measurement package according to (4.6) can be approximated with 
more degrees of freedom compared to the PCA. Therefore, for each given 
point cloud 


P[n] = |p,[n],---,Pm[n],---,Pain]] E RM, (4.37) 


the task is to find a 2D curve f(r) by linear combination of F B-splines, 
which fits the given point cloud best in a least-squares (LS) sense. For 
simpler notation, the dependency on the discrete time step n will be 
omitted for the remainder of this subsection. It is only important to keep 
in mind that the presented approximation must be applied for each point 
cloud of the reference signals P,[n] and delay signals P [n]. 

The fitting problem of the spline curve to the M given 2D points can 
be reformulated in the matrix form 


i J 1" e o] | Cx,0 Cy 0 
: : = : = : : : +N, 
1m PaM boi (Tas) a brail m) Ex,F-1 Ĉy,F-1 
PTeRM~x2 B(r)ERM*xF CeRF*? 
(4.38) 


with the approximation error N € R™*? and the corresponding support 
points 


r= [rye ru, rm €{R|0<r,, <1}. (4.39) 


The problem (4.38) can usually be solved using an LS estimator. However, 
to build the Moore-Penrose inverse of B(r), besides a set of B-splines, a 
corresponding support point r,,, for each point within the point cloud is 
necessary. In simplified terms, this means that the support points r need 
to be estimated first. Unfortunately, the determination of these support 
points influences the precision of the approximation and the complexity 
of the curves x(r), y(r) a lot. In order to get smooth curves with the 
least possible curvature, the distance between two support points should 
be proportional to the distance of the corresponding points along the 
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spline curve, which in turn is only available after the support points are 
determined. 

Consequently, the problem of finding the support points and the spline 
curve is solved jointly by the iterative method of Wang et al. [127]: 
0). 


rí initial support points 


CHD = (BT(r)B(r®)) T BT(r®). PT 4.40 
(4.40) 


pth) — arg min |[B(r) Ce) _ prp 


Finding the initial support points r is realized by pre-ordering the 
points p„ using a projection onto an ordering curve. Furferi etal. [37] 
presented an algorithm for ordering the points by piecewise progressing 
along the line. The objective is to find ordered contour points along the 
point cloud. In this approach, starting from a random point of the cloud, 
the direction in which to search the next contour point is determined 
using a PCA of the points in proximity. The radius of the circle to differ- 
entiate between points lying in proximity and further distanced points is 
gradually increased until the eigenvalues of the PCA reach a predefined 
ratio. After the radius is fixed, the next contour point is calculated to be 
the intersection of the first principal component with the circle. While the 
idea of determining contour points to get an ordering curve was taken 
up in this work, the PCA-based algorithm to find them could not keep 
the robustness requirements. 

Therefore, a simpler but more robust approach, which selects contour 
points from the cloud, is developed. In the following context, contour 
points p,,, are distinguished from original points by a tilde. Firstly, the 
number of contour points has to be set. Examination of different point 
clouds from the simulated data sets has proven that it is favorable to 
select ten contour points, thus, this hyperparameter is fixed to ten. Since 
the trajectories do not have loops and are only slightly curved even for 
higher signal dynamics (see Fig. 4.2), the two points with the maximal 
distance between each other represent the start and end points 


1<i j <M. (4.41) 


{B1 Pio} = arg Max |p: = P,| 27 
Pi>Pj 
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Figure 4.7 Visual example of iterative 2D B-spline approximation. 


The remaining eight support points p,,, are selected in a way that they are 
distributed as uniformly as possible over the maximum point distance: 


z . z m-— 1 
P, = arg min Ip; — Bill, — —g max , 2<m<9, (4.42) 
Pj 


with 
dmax = Bio = Pill, : (4.43) 


After all contour points are found, their support points ?, are esti- 
mated using the normalized cumulative Euclidean distance 


7 1 K MEN 444) 
r = : = . 
= BB, Ejea lP; — Ball, mel 


The scaling of the support points ? „ can be arbitrarily chosen because the 
points to approximate are only given as 2D vectors, but the scaling needs 
to be adapted to the interval spanned by the knot vector. For simplicity, 
the scaling of the support points and correspondingly the knot vector 
are set to the interval [0, 1]. Subsequently, a B-spline approximation of 


73 


4 Novel Separation Methods Based on Signal Dynamics 


the contour points with the support points 7,,, yields the ordering curve 
PO (r), which in turn is employed to find the initial support points r ® 
by orthogonal projection of the 2D points onto the ordering curve. Due to 
the analytic derivatives and formulations of the B-splines, this projection 
can be solved analytically by intersection of the normal vectors with the 
points. 

Figure 4.7 depicts the results of finding the contour points p,,,, the 
ordering curve, and the final spline curve. On the right side the con- 
vergence behavior of the summed squared orthogonal distances can be 
observed. Due to well-chosen initial support points, the spline curve 
converges after only a few iteration steps. The residual mean squared 
error (MSE) should not be significantly less than the standard deviation 
of the measurement noise, as this would indicate overfitting. In this case, 
the degrees of freedom F have to be reduced by using a knot vector with 
fewer elements. For both the ordering curve and the iteratively approxi- 
mated spline curve, the B-spline order and the knot vector are set to | = 4 
and r, = [0, 0.5, 1]", respectively. 


4.3.3 Interfering signals estimation 


From the discussion about the misalignments of the point clouds in 
Section 4.1.3 it is clear that determining these misalignments yields the 
differential interfering signals in its vector form p^ [n]. However, there 
are several problems to be solved. Firstly, the overlap between two point 
clouds p[n], p[n+1] is influenced by the interval t, which in turn depends 
on the interval of the time delays. Hence, the length of the overlap is 
not known beforehand. Secondly, the misalignment cannot be calculated 
from the point clouds directly, because the distribution of the sampled 
time delays in the interval t may be non-uniform or even contain gaps 
leading to points p„ [n] with no corresponding point p, [n + 1] in the 
subsequent point cloud. This disqualifies the application of point set 
registration approaches such as iterative closest point since they require 
existing direct point correspondences. To circumvent this issue and to 
be more robust against AWGN, the shapes of the points clouds are first 
approximated by using B-splines as explained in the previous subsection, 
instead of calculating the misalignment between the point clouds. For 
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each time step n, this results in two trajectories represented by the spline 
curves p,,(r) and p,,, ,(r). As the shapes only overlap partially, finding 
the misalignment results in the problem called partial shape matching (see 
[141] for further reference). 

To solve this problem, the trajectories p,,(r), P,,,.,(”) are densely sam- 
pled and the misalignment is determined by minimization of a mean 
squared distance-based quality measure. Because the curves are only 
partially overlapped, the crux in the quality measure is to classify which 
points have a corresponding point in the other trajectory. To this end, 
the advantage that sampling the B-splines creates ordered points is ex- 
ploited. Considering the direction of the trajectory, the end point of p,, (r) 
which is guaranteed to have a correspondence on p,, , , (7) is easy to de- 
termine. Following from this approach, the misalignments used in the 
optimization are reduced to the differences 


BS, (r) = Pn41 (r1) — Pn (0) (4.45) 


at the sampled support points r,. Without loss of generality, Equation 
(4.45) assumes that the end point with guaranteed correspondence is 
p,,(0), otherwise the control points of the spline curve can simply be 
reversed. Thus, the 2D optimization is simplified to a 1D search. 

Given N densely sampled support points r, € [0,1] and the misalign- 
ments from (4.45), the final estimation for the differential interfering 
signals can be calculated via 


N 
BA In] = argmin min |p, 10) — Bary) — BEM), 
PF, (Ti) i=i 


(4.46) 


Equation (4.46) minimizes a quality measure formed by accumulating 
the Euclidean distances of the respective nearest neighboring points, 
where the starting index i is responsible for classifying which points 
have a correspondence. 

The last step for estimating the interfering signals is the numerical 
integration of p^ [n] followed by offset correction. To this end, the rect- 
angle rule is applied and the resulting cumulative sum is searched for 
local maxima. Since the interfering signals are bandpass signals, the 
average calculated over one period should be zero. This is exploited by 
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Figure 4.8 Example estimation of the interfering signals. Upper left: all point clouds. 
Upper right: B-Spline approximation of a single point cloud. Lower left: partial shape 
matching between two trajectories. Lower right: estimated differential interfering signals 
sA(t,) and their offset corrected integration. Figure adapted from previous publication 
[141]. 


calculating and subtracting the average value of the range between the 
found local maxima. 

An example of the entire algorithm consisting of the point cloud gen- 
eration, the B-spline approximation, the partial shape matching, and the 
numerical integration with offset correction is shown in Figure 4.8. For 
this illustration, a measurement package was extracted from the data 
set D.i,,1- After clipping the time domain containing the target signals, 
the point clouds of the reference signals are depicted for each time step 
n in the clipped time domain (upper left). Although not shown, each 
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point cloud is subsequently approximated by a spline curve (upper right), 
which in turn are used to calculate the misalignment for each time step 
n (lower left). On the lower right, the estimated differential interfering 
signals are plotted in blue and their numerical integration in red. The 
shaded area marks the range, where the average for the offset correction 
was calculated. Due to the solution in the 2D area, the interfering signals 
are estimated in form of analytic signals, but these can be easily converted 
back to the real-valued interfering signals by simply taking the real part. 

The last point to mention is which signals are to be used for the point 
cloud generation since the B-spline-based approach works with both the 
reference or delay signals. There are two options to deal with it: either 
the interfering signals are estimated once using the reference signals and 
once using the delay signals followed by averaging, or the point clouds 
are generated by concatenating the reference signals and the delay signals 
leading to P[n] = [P,[n], Pg{n]]. For easier notation, the hyperparameter 
set 


OPS = {Ar,PCgen}, PCgen € {individual, concatenated} (4.47) 


is introduced, with the categorical variable PCgen which decides whether 
the point clouds for the interfering signals estimation should be concate- 
nated or the point clouds of the reference and delay signals are generated 
and processed individually followed by an averaging. 


4.3.4 Joint shape fit and interfering signals estimation 


The B-spline-based SDM proposed in the previous two subsections has 
two drawbacks. Firstly, using an individual B-spline approximation for 
every time step t,, yields slight deviations in the shapes of the trajectories 
due to the noise n leading to further errors when solving the partial 
shape matching problem. Secondly, the solution of the partial shape 
matching requires a lot of computational effort, since the approach is 
based on sampling the splines, where the high precision requirement of 
the misalignment requires a high sampling density. Therefore, a novel ap- 
proximation method, named Joint B-Spline approximation, is presented 
that considers the knowledge that the shapes of the trajectories at dif- 
ferent time steps are constant except for the misalignment. By using a 
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joint shape fit and misalignment estimation, a global optimum can be 
reached. To accomplish this, the fitting problem (4.38) is extended by the 
misalignment, leading to the formulation 


P"[n] B(r[n]) On C N[n] 
= . , (448 
Paar Br tint bäi] Im rl P 
En a gr in I RE 
PT[n]er?Mx2 B(FJER2MX(FHN) GcRF+Vx2 


with the ones vector 1,71 € {1}”7, the zeros vector O7, € {0} and the 
transposed misalignment p4[n]. Instead of fitting the B-splines to only 
the point cloud P[n], the same spline curve needs to approximate both 
consecutive point clouds P[n], P[n + 1], with the second point cloud 
being shifted by the estimated misalignment. Inserting the modified 
matrices P[n], B(#) and © into the iterative approach (4.40) results in 
both the final spline curve and the differential interfering signals at the 
time step n. Finally, the numerical integration with offset correction as 
already explained in the previous subsection is repeated to obtain an 
estimate for the interfering signals. 

However, the quality of the initial support point estimation deteriorates 
with increasing offset and the convergence is not guaranteed anymore 
as Figure 4.9 illustrates. Since the convergence rate is rather low and the 
LS estimation in (4.40) can be interpreted as the direction of a steepest- 
descent approach, the learning rate could be increased. To this end, a 
golden-section search was implemented, but the results showed that 
manipulating the learning rate during the iterative optimization makes 
the convergence unstable. For this reason, the line search was omitted 
when applying the joint B-spline. Nevertheless, good results could be 
achieved by increasing the iterations by factor 10 compared to the B-spline 
approximation of each individual point cloud. Overall, the computational 
effort was still below that of the individual approach, because the effort 
for the partial shape matching could be saved. 


4.3.5 B-spline-based direct time delay estimation 


Equivalent to the direct TDE using the angle of the first principal compo- 
nent vector, the spline curves can also provide a tangent line. Furthermore, 
due to the estimation as a curved line, it is possible to determine distinct 
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Figure 4.9 Visual example of iterative joint B-spline and misalignment approximation. 


angles for each measurement within the package, instead of only one 
orientation per point cloud. Because the B-splines can be derivated ana- 
lytically, the angle calculation is fairly simple. The following steps must 
each be applied once to the point clouds resulting from the reference 
signals and the delay signals. At first, for each point p„ [n] the support 
point with the least Euclidean distance 


rIm,n] = arg min [[p,, (r) — Pmlnlll, (4.49) 
re[0,1] 


has to be computed, with p, (r) being the spline curve approximated to 
the respective point cloud P [n]. Secondly, the analytical derivative 


dp,, (7) 
vgs[m, n] = es (4.50) 


dr r=r([m,n] 
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is calculated at the support points obtained in the previous step. As the 
spline curve is a weighted sum of the B-splines, the analytical derivative 
becomes 


dp,, (7) = [6 | s db; (r) 


g : 451 
i Cy,i dr l ) 
where the derivation of the B-splines needs to be substituted by the 
expression given in (4.36). Finally, the angles 


gsm, n] = atan2 (e} vasim, n),eTvps[m, n]) (4.52) 


of both the reference and delay signals are sorted per measurement 
index m, resulting in drm [n] and bam [n]. From these angles, again using 
the same procedure already applied for the PCA-based direct TDE in 
Subsection 4.2.2, the time-delay difference 7,,, is determined for each 
measurement index m. 

Since this direct TDE is only based on the spline curves, it can be used 
for both the individual B-spline and the joint B-spline approximation. 
Furthermore, similar to the B-spline-based interfering signals estima- 
tion, the hyperparameters specify the signal dynamics and whether the 
point clouds for reference and delay signals should be concatenated or 
processed individually (cf. the set defined by (4.47)). 


4.3.6 Discussion 


The B-spline-based point cloud processing has several advantages but 
also disadvantages to the PCA-based SDMs. The first point to mention 
is that the B-spline-based SDM can also be applied if only either the 
reference or the delay signals are available. Furthermore, due to the 
estimation via the differential signal, which needs to be integrated before 
the interfering signals are acquired, the B-spline-based SDMs have a 
better robustness against noise. The next advantage lies in the dynamic 
behavior if the time delays are estimated directly from the tangent lines 
of the spline curves. Here, the B-spline-based point cloud approximation 
allows us to obtain an estimation for each signal within the measurement 
package rather than just an averaged value per package. The last point to 
mention is that the requirements on the point clouds are less restrictive. 
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Figure 4.10 Results of TDEs using the B-spline-based SDM for the simulated data sets 
Dsim,1> ey Dgim,4- 


The variations of the time delay within a package may be larger on the 
one hand and on the other hand the distribution of the points within the 
point clouds are not required to be uniform. 

However, the computational effort is considerably higher than for 
the PCA-based processing and due to the more complex algorithm, the 
method is more susceptible to outlier signals, e.g., zero signals due to 
loose contacts, fluctuating time shifts due to wrong trigger signals, etc. Es- 
pecially, the performance of the joint-B-spline approach depends on the 
quality of the signals, since the convergence is not guaranteed. Therefore, 
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outliers would have to be filtered out beforehand to get reliable estimates. 
The last expected disadvantage is the robustness against model errors. 
Because the algorithm is based on the invariance of the shape against 
time shifts of the target signals, any shape manipulating effect, such as 
varying interfering signals, should have a strong impact on the misalign- 
ment estimation. Of course, a changing amplitude of the target signals 
could also have a negative effect, but during the measurement system 
identification this effect was not observed and is therefore considered as 
unlikely. 

In order to investigate the discussed properties, the joint-B-spline- and 
B-spline-based SDMs are applied to the four simulated data sets leading 
to the results plotted in Figure 4.10. Same as for the PCA-based approach, 
the signals are separated into packages with Ar = 50 ns, but due to the 
higher computational effort, the overlap was reduced to allow the calcu- 
lations in reasonable time. From this follows the staircase-like curve that 
can be observed in the upper and lower left of Fig. 4.10. Furthermore, 
the estimate of the curvature of the splines at the endpoints of the trajec- 
tories deteriorates, resulting in poor B-spline- and joint-B-spline-based 
direct TDE at these points, which explains the discontinuities in the TDE 
for signals at the edges of the signal packages. Regardless, the figure 
shows that the spline-based approaches work in general but, as already 
suspected, show stronger deviations than the PCA-based approaches in 
case of model errors. The joint-B-spline-based processing can cope better 
with model errors, but still not as good as the PCA-based processing. In 
return, however, it is also allowed that strongly varying time delays are 
present in the measurement package for the direct TDE (see the upper 
and lower right in Fig. 4.10). Since the simulated data sets are generated 
with a high SNR, the robustness against noise cannot be investigated 
in this evaluation. Therefore, this point will be discussed later in the 
experimental results. 


4.4 Global interfering signals compensation 
In practice, it would not be usual to recalculate the interfering signals 


for each measurement package without taking into account the informa- 
tion given by already existing estimates of the interfering signals. While, 
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due to the interfering signals being slightly temperature dependent, it 
is beneficial to have an individual estimation per package, the robust- 
ness can be further improved using a globally valid estimation for all 
measurement packages. To this end, all estimations obtained from evalu- 
ation of the different packages are averaged after outliers were filtered 
out. Subsequently, all measurement signals can be compensated globally 
using the thus obtained interfering signals estimation, followed by the 
conventional zero-crossing-based TDE. This approach is referred to as 
global compensation SDM later in the results chapter, whereas the com- 
pensation of each measurement package with its individual estimation 
will be called local compensation SDM. 

As already mentioned in the discussion of the B-spline-based process- 
ing, outlier points in the point clouds can greatly impact the quality of the 
estimation. To remove bad estimations, a new histogram-based outlier 
detection algorithm is introduced in the following. Given M different 
estimated interfering signals s;[n, m], in the first step, a histogram 
M & a nn 

h ee (453) 
0 siin, Mm] é fb; ,b;,,) 


4? “44 


is calculated for each discrete time index n, where h? denotes the number 
of samples that fall within the interval [b, Bag These intervals, also 
called bins, are determined for each time index n separately. For this, 
the next lower power of ten to Scott’s normal reference rule [106] is 
chosen as bin width. Then, the intervals are uniformly distributed over 
the maximum interval | min,,, 5,[n, m] , max, s,[n, m]], determining the 
time-index-specific interval limits b;’ and thereby the number of bins. 
In other words, the interfering signals sampled at the time index n are 
treated as a collection of values, which is then used as input for an 
automatic binning algorithm. Subsequently, the bins are applied to get a 
histogram of the collection of values. As this procedure is repeated for 
all time steps, a set of histograms is eventually obtained. With the help 
of these histograms, a function 

he 8 € [by bral 
AP:ERHN, sr LS: (4.54) 


ae ae] 
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can be defined for each time step, which assigns to a continuous value 
s the number of observations that lie in the interval in which the value 
also lies. 

In the second step of the outlier filtering, the histogram mapping 
function Å” is applied to create a quality measure that indicates which 
signals are less likely to be outliers and which are more likely to be 
outliers. Since a signal sample being in the bin with the highest number 
of observations indicates that it belongs to the majority, a higher value 
means a signal that is less likely to be an outlier. Obviously, the signals 
do not consist of only one sample. Therefore, the number of observations 
for each signal 5, [n, m] are accumulated over the samples in the quality 
measure 


N 
Jin] = Sh" (sin, M)). (4.55) 


Finally, a threshold determines for which J[m] the signals are con- 
sidered outliers. An automatic threshold estimation can be found by 
sorting J|m]. If the first and the last point of the found quality measure 
are interpreted as two 2D points connected by a straight line, the point on 
the quality curve with the largest Euclidean distance to this straight line 
is taken as the threshold, i.e., the strongest bend in the curve is used as 
separation point. The last task is to remove the outliers and average the 
remaining signals to get a globally determined estimate of the interfering 
signals 5,(t,, ). The results of using this globally determined estimate will 
be shown in the experimental results chapter. 
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A completely different approach, which is not based on compensating 
the interfering signals but on estimating the time delays robustly against 
interfering signals, is presented in the following. In principle, features 
are extracted from a time-frequency representation of the signals, from 
which the time-delay difference is subsequently estimated with the aid of 
a regression. To present this method, the fundamentals of the used time- 
frequency representation are given in the first section, followed by an 
overview of the algorithm in the second section, the feature design in the 
third section, and a feature subset selection in the fourth section. Finally, 
the last section of this chapter examines different regression models and 
the corresponding training approaches for both the supervised and the 
unsupervised case. 


5.1 Fundamentals on analytic wavelet packets 


The analytic wavelet packet transform (AWPT) belongs to the class of 
multirate filter banks, just like, e.g., the discrete wavelet transform (DWT) 
[83]. The applications of these multirate filter banks include, but are 
not limited to, speech processing, image processing, data compression, 
and denoising as they are fast to calculate and invertible. Most of these 
applications take advantage of the fact that the multirate filter banks 
create a time-frequency representation of the input signal that provides 
information about the signal energy in a certain time and frequency win- 
dow. Since the time-frequency resolution of the DWT is fixed to have a 
high frequency resolution at low frequencies and a high time resolution 
at high frequencies, the more flexible wavelet packet transform [19] was 
introduced. However, the drawback of the real-valued multirate filter 
banks is the missing shift-invariance due to the downsampling operation 
in the algorithm. Although many variants retaining the shift-invariance 
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were investigated in the early 2000s, such as the undecimated wavelet 
transform [116] or the double density DWT [108], this thesis uses the 
analytic wavelet packet transform, which was published independently 
by Selesnick [110] and Weickert et al. [128] at about the same time. There 
are two main reasons for this decision. Unlike the undecimated wavelet 
transform, the amount of data does not increase at each level, and more- 
over, the use of analytic basis functions allows the determination of the 
time shift based on the phase of the complex-valued coefficients just as 
in the Fourier transform. 

The following subsections explain the structure of the AWPT and the 
principles to design the necessary tree-like multirate filters and filter 
impulse responses. The presented fundamentals are summarized from a 
collection of publications on the subject [24, 58, 107, 109-111, 128, 144]. 


5.1.1 Analytic wavelet packets in the context of 
time-frequency representations 


Basically, the AWPT can be classified as a basis transform for discrete sig- 
nals. Representing signals in other basis systems to extract information 
is a common technique. The most prominent basis transform is probably 
the Fourier transform which allows the analysis of the frequency com- 
position of signals. In the discrete version of the transform, harmonic 
oscillations with different frequencies form an orthogonal basis, i.e., the 
inner product, which measures the similarity of two functions, of two 
distinct basis functions is zero. Other known basis transforms are the 
wavelet transform or the STFT, whose main difference to the Fourier 
transform is the compact support of the basis functions in the time range 
allowing a time-dependent statement about the frequency composition. 
Whereas the STFT can be interpreted as the Fourier transform of a time- 
windowed signal, the wavelet transform is defined as the inner product 
between a signal and a time-stretched and time-shifted mother wavelet. 
In both cases, the form of the basic functions is determined by the window 
function and the wavelet, respectively. Even though the time-frequency 
resolutions can be adjusted by the window length for the STFT or the 
bandwidth or center frequency of the mother wavelet, the general traits 
regarding the time-frequency resolution of these transforms—the STFT 


86 


5.1 Fundamentals on analytic wavelet packets 


has a constant time-frequency resolution and the wavelet transform has 
a better frequency resolution at lower frequencies—are not changeable. 
Furthermore, the product of the frequency resolution coș and the time 
resolution g, is larger or equal to the Gabor limit 1/(47), where equality 
only holds for the Gaussian window function. 

In this context, the DWT is the implementation of the real-valued 
wavelet transform for discrete signals by using filter banks with different 
sampling levels to replicate the time stretching of the mother wavelet. 
Therefore, it shares the same traits—better frequency resolution at lower 
frequencies and oscillating coefficients in some cases even for signals with 
a constant frequency. An extension of the DWT to use complex-valued 
basis functions equivalent to the DFT is called the dual-tree complex 
wavelet transform. This approach is already very close to the AWPT. The 
only extension missing to get the AWPT is the adaptive time-frequency 
resolution, i.e., different time-frequency resolutions can be set arbitrarily 
for different frequency levels within the Gabor limit. Basically, however, 
all presented time-frequency representations are still only basis trans- 
forms where the basis functions have certain properties, e.g., bandwidth, 
center frequency, time compactness. 

Mathematically, a basis transform of a discrete signal s[n], e.g., the DFT, 
can be formulated as finding the coefficients c[k] needed to decompose 
it into the sum 

K-1 


sin] = I clk] yin]. 6.1) 


k=0 


Although the signal s[n] may in principle be real-valued or complex- 
valued, this thesis only considers real-valued signals due to the origin 
of the signals from ultrasonic measurements. To get the coefficients, the 
Gram matrix of the new basis system 


® = Basis {y [n] , k=0...K—1} (5.2) 


is inverted and multiplied by the vector resulting from building the in- 
ner product of the signal with each basis function x, [n]. Note that the 
decomposition in (5.1) is only unique if the number of basis functions K 
is equal to the number of discrete samples N. For K > N the function 
system, called frame or dictionary in this case, is overcomplete and the 
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coefficients are determined by imposing additional constraints. Other- 
wise (K < N), the system only spans a subspace of all possible signals 
leading to the situation that only an approximation of the signal can be 
found. In all cases, the inverse transform is simply building the sum (5.1). 

A special case is given when the basis is orthonormal. Since the inner 
product of distinct basis functions equals zero and the signal energy of 
the basis functions equals one, the Gram matrix becomes an identity 
matrix. Therefore, the coefficients can be determined directly through 
projection of the signal onto the respective basis function 


N-1 
clk] = (s[n), d.[r]) , = > sir] vy ln. (5.3) 
n=0 


Even if the basis is only orthogonal, normalizing the signal energy of 
each basis function creates an orthonormal basis. Furthermore, for or- 
thonormal bases applies: given the basis function Y, [k] the coefficient 
c[k] not only provides information about how much signal energy is in 
the time window around the mean time es but also at the same time, 
due to Parseval’s theorem, how much signal energy is around the mean 
frequency f p- Here, the mean time t g, the mean frequency f g the time 
duration o, p, and the bandwidth o; ;, are defined by the relationships 


= = a a = - , Vln]? 
A DE N = I (nts — t, ———, (54) 
n=0 
F N ong a = (= u ) [T in]? 
k” a Ope ; , 
np=0 mli l ng=0 Nt, EZAZ 


(5.5) 


where Y, [n] denotes the DFT of the basis function y, [n]. 

Summarizing the above-mentioned properties of basis transforms, 
where the objective is to get information about the time-frequency energy 
distribution of signals, the AWPT is designed to have the following 
properties: 


= The transform should use complex-valued analytic basis functions 
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vein] = vIn] —j yin], with yin] = Hl [n]}, 6.6) 


where the real and the negative imaginary part build a Hilbert 
pair, i.e., the negative imaginary part is the Hilbert transform of 
the real part. Due to this property, the shift-invariance of the trans- 
form is preserved, although the downsampling operations are still 
included. 


In order to calculate the coefficients efficiently as inner products— 
later it will be shown that the AWPT applies discrete convolutions 
to perform this task—the basis functions wr [n] should be con- 
straint orthonormal, i.e. 


(yr in], oF in}, = 056,5, 
(vy), = 058 0SE7SN-1, 67 
(yenyen) ~o, OSkSN-1, 


Note that if condition (5.6) does perfectly hold, the inner product 
between the real and imaginary part is zero. However, in prac- 
tice, as it will be shown later, the filter design using finite impulse 
response filters can only approximately design the real and imagi- 
nary parts as a Hilbert pair, since the Hilbert transform requires 
an infinite impulse response. As the set of real-part and imaginary- 
part basis functions each form a complete basis system, the inner 
product between real and imaginary parts from different frequency 
bands k is not required to be zero. 


In order to perform filtering or other processing tasks, the AWPT 
needs to have a numerically stable inverse transform. 


Each signal has its own time-frequency characteristic. Therefore, 
a rigid time-frequency resolution is not optimal for every signal. 
To this end, the time-frequency resolution of the AWPT should be 
signal adaptive, where the resolution is chosen to be best under 
certain boundary conditions such as minimal entropy [19]. 
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5.1.2 Dual-tree approach for multirate filter banks 


A multirate filter bank such as the DWT or wavelet packet transform 
is built as a full binary tree, i.e., each node can have exactly two or 
zero children, where a node without children is called a leaf. In this 
context, a node consists of a lowpass or bandpass filter followed by a 
downsampling by two. Due to this, successive filtering steps lead to a 
dyadic reduction of the sampling rate, which is why these structures are 
also called multirate filter banks. Furthermore, the filtering influences 
the frequency windows that the coefficients represent, i.e., the frequency 
characteristic of the basis function corresponding to anode is determined 
by the sequence of filters that leads to the respective node. In principle, 
two main properties of this frequency characteristic are set. Firstly, as 
each node reduces the sampling rate, the time resolution is halved while 
the frequency resolution is doubled for each additional filter pair that 
is passed through. Secondly, the order of the passed filters determines 
the center frequency, e.g., a bandpass followed by a lowpass yields the 
center frequency fy = 5/8fn, with fy denoting the Nyquist frequency. 

In summary, given the input coefficients c[n] of a node, the coefficients 
of the subsequent lowpass and bandpass children are calculated by 


crpln] = (er). ,,- (58a) 


cpgp[n] = (c[fi] u Ip |R]) (5.8b) 


A=an, 
with the impulse responses of the lowpass filter g,,p[n] and the bandpass 
filter ggp[n]. Since the filtering results in two sequences of coefficients, 
each of which is subsequently downsampled by two, the number of 
coefficients at each level remains constant. Note that the information 
lost due to discarding half of the samples after each filtering can only be 
recovered if the impulse responses are designed to hold certain properties, 
which will be discussed in the next subsection. 

An example of a full multirate filter bank representing the AWPT is 
depicted in Fig. 5.1. Additionally to the already mentioned relationships, 
a dual tree responsible for calculating the imaginary part of the coeffi- 
cients is included. This approach can be explained as follows. Recursive 
application of (5.8) leads to the coefficient sequences of the nodes, which 
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Figure 5.1 Dual-tree representation of the analytic wavelet packets. The filter pairs are 
rearranged to obtain analytic basis functions and to order the center frequencies ascending 
with the subband (see Section 5.1.3). 


represent the inner products of the corresponding basis functions with 
the input signal. However, since the calculations are performed all in the 
real-valued domain, using a single tree yields only the real part of the 
coefficients c*°[k], or in other words, only the basis functions pie [n] are 
implemented by a single tree. The dual-tree approach adds a second tree 
and slightly modifies the impulse responses to guarantee that the implicit 
basis functions of the original tree ape [n] and the dual tree yin [k] build 
a Hilbert pair—the order of the filter impulse responses is modified to 
get the correct behavior with respect to the center frequency and analytic 
basis functions, as explained in the following subsection. With this, the 
final complex-valued coefficients can be obtained by 


c[k] = (s[n], V(r] - j vp [n]a 
= (s[r], YE In), +j (sir), VET) 
= eek] + jcl™[k] . (5.9) 
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Figure 5.2 Time-frequency resolutions of two different tree structures. The standard 
time-frequency resolution of the DWT is shown on the right side. 


Another important aspect that has been neglected so far is the nomen- 
clature of the coefficient indices. Until now the coefficients were only 
indexed by the discrete frequency step k. However, a detailed examina- 
tion of the multirate filter bank shows that each node is uniquely defined 
by its level and subband. Furthermore, each node contains multiple 
coefficients due to the time index n. 

Firstly, the node nomenclature used in this thesis—and in the previous 
publication [144]—is explained. Each node is described by the tuple 


(s,,k,) EN?, with s, >0,k, € 0,2-1], (5.10) 


where s,, denotes the level of the node and k,,, denotes the subband within 
this level. Starting from the root with the tuple (0, 0), the following nodes 
are determined recursively by a parent-child relationship, in which the 
children are defined as 


lowpass child (s„ + 1, 2k„), (5.11a) 
bandpass child (s„ + 1,2k„ + 1), (5.11b) 


if the parent node is given by (s,,, k,,). 

Secondly, the time index n of the node coefficients is set in relation to 
the time-frequency window that corresponds to each individual coeffi- 
cient. While the mean frequency is completely determined by the node 


identifier (s,,, k„), the mean time is determined by the time index of the 
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coefficient sequence. Unfortunately, the coefficients cannot be arranged in 
the form of a matrix, because the coefficient sequences c, p [n] from dif- 
ferent nodes are not always of the same length. Therefore, the coefficients 
are vectorized into a linear arrangement by simple concatenation of the 
different coefficient sequences and indexed by k, which in this context 
specifies both the mean frequency and mean time. Due to this, a map 
needs to be stored containing the level, the subband, and the time index 
for each index k. Two examples for different tree structures, their coeffi- 
cient notations c, | ;,[n], and their time-frequency resolutions indicated 
by rectangular windows are depicted in Fig. 5.2. Due to the limited space 
the time duration is shortened to include only eight samples, which can 
be concluded from the number of coefficients in the first level (s,, = 1) 
being four. 


5.1.3 Filter design and arrangement 


The desired properties of the AWPT stated in Subsection 5.1.1 and the 
existence of its inverse transform are strongly dependent on the design 
of the filters and their arrangement. Since the focus of this thesis does not 
lie on the filter design, only the fundamental approaches are presented. 
For further reference see the works of Daubechies [24], Selesnick [107, 
109] and Kingsbury [58]. 

In the case of the conventional DWT without a dual tree and with 
non-analytic basis functions, Smith and Barnwell [115] have proven that 
a bandpass-lowpass filter pair with the relationship 


Ipplr] = (—1)"gpp[N — n] (5.12) 


allows a perfect reconstruction, which is in this context the inverse trans- 
form, if the perfect reconstruction condition 


Gip(z)Gyp(z *) + Gpp(—z ")Gp(—z) = 2 


l (5.13) 


S sul sreli + 2n] = fr] = $ um 


is fulfilled. Such a filter pair is called conjugate-quadrature-filter and 
leads to orthogonal filters. The reconstruction or synthesis filters h[n], 
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which are used in the mirrored synthesis tree not shown in this thesis, 
can then be derived from the lowpass filter according to: 


hipin] = grp N- n], (5.14a) 
h 1)” ttg pin]. (5.14b) 


In summary, a set of four filters is necessary to perform the DWT and its 
inverse transform. Note that the perfect reconstruction condition (5.13) 
still contains degrees of freedom. One option to design such a lowpass 
filter is using Daubechies’ approach 


Grp(z) = (1+ 2)"Q(2), (5.15) 


leading to the so-called Daubechies-N wavelet, where N denotes the 
resulting length of the filters [24]. In this context, the number of vanishing 
moments L has to be set and a polynomial Q(z) has to be designed, which 
can be done by a spectral factorization. Furthermore, there are several 
more design methods, such as Symlets with improved symmetry or 
Coiflets where the resulting scaling function also has vanishing moments 
[82]. Whereas the standard Daubechies-N wavelets lead to a minimum- 
phase system, the Symlets and Coiflets are nearly linear phase systems. 

Until now, only the filter design for the real-valued DWT has been 
addressed. However, the extension to the AWPT requires further filters 
for the second tree and its reconstruction tree. Additionally to the perfect 
reconstruction conditions, the resulting wavelets need to build a Hilbert 
pair. As proven by Selesnick, this can be realized if the filters of both 
trees satisfy the half-sample delay condition 


gpl] = gppln — 0.5]. (5.16) 


Since the filter impulse responses are discrete, this condition is not as 
trivial as it may seem and can only be solved approximately. Three promi- 
nent design approaches are the q-shift wavelets [58] and the linear-phase 
biorthogonal wavelets [59, 60] introduced by Kingsbury as well as the 
common factor solution presented by Selesnick [109]. For the latter, the 
vanishing moments of the resulting wavelets and the approximation 
accuracy of the half-sample delay condition can be set independently by 
two parameters. 
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Figure 5.3 Reason for exchanging the bandpass and lowpass filters after bandpass node 
in multirate filter bank (adapted from [56]). 


However, Selesnick et al. [111] realized that even if the half-sample 
delay condition (5.16) is met perfectly, the resulting wavelets are only 
analytic if enough filters are passed through. This problem can be solved 
by using different filters in the first level, which satisfy 


Oui [n] = dr [n- 1]. (5.17) 
Here, the different filters are highlighted by a subscripted 0. Because this 
condition is a lot easier to satisfy, ordinary Daubechies wavelets can be 
used in the first stage of the real-part tree and their one-sample delayed 
version in the imaginary tree. The corresponding filters for the synthesis 
trees can be calculated again by the relationships (5.12) and (5.14). In 
summary, a total of eight analysis filters and eight synthesis filters is 
necessary to implement the AWPT and its inverse. 

Finally, the filter arrangement has to meet a few criteria to get analytic 
basis functions and to order the center frequencies f, ascending with 
the subband k,,. As Figure 5.3 illustrates, the downsampling shifts the 
bandpass parts sgp[n] into the lower frequency band, but with swapped 
lowpass (L) and bandpass (B) components. To account for this effect, the 
bandpass and lowpass filters must be swapped in the next stage (see 
Fig. 5.1). After another bandpass, the filters are swapped again. Due to 
the same effect, the filters of the imaginary- and real-part tree have to 
be swapped after the first bandpass (see Fig. 5.1). Otherwise, the basis 
functions would have only energy at negative frequencies. The last point 
to mention is that after reaching a node in the dual tree, where the basis 
functions are already analytic signals, the subsequent nodes in both the 
real and the imaginary trees need to have the same filters [110, 128]. For 
simplicity, the filters of the real-part tree are used in this case. 
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5.2 Algorithm design 


The proposed approach for the robust time-delay difference estimation, 
in the following called the feature-based regression method (FRM), is 
inspired by classical machine learning methods. To this end, the data 
set needs to be split into a training and test set. Subsequently, feature 
vectors are extracted and, after a reduction of the dimensionality, used 
for a regression to estimate the time-delay differences. 

There are two general classes of this approach: supervised and unsu- 
pervised. In both cases, the objective is to learn from a given training 
set how to estimate the correct time-delay difference given a measure- 
ment signal each for the reference and delay signal. However, in the 
supervised case, the corresponding ground truth, also known as labels, 
for the time-delay difference of the signals in the training set must be 
known. In contrast, unsupervised machine learning methods only need 
the training set without any corresponding labels. In practice, this means 
that a calibration in a scenario with known ground truth is necessary for 
the supervised approaches, whereas the unsupervised approaches can 
collect the data during operation and learn the regression from it. As a 
result, unsupervised approaches are closer to the SDMs or to conven- 
tional TDE methods and have higher acceptance in practical applications. 
For this reason, even though both supervised and unsupervised versions 
of the FRM are proposed, the focus lies on the unsupervised version. A 
supervised version of this class of approaches is already published in 
an own previous publication [145]. This chapter takes up the approach 
contained in [145] and extends it by adding further selection methods, 
regressions, and a solution approach for the unsupervised training. 

The principle of the supervised FRM is outlined in Fig. 5.4. By omitting 
the label connections to the training of the feature selection and the 
regression model, the approach is modified to get the outline of the 
unsupervised FRM. Since the proposed feature extraction algorithm is 
based on the AWPT, the resulting number of features is much too large to 
get a stable regression. Therefore, a subsequent feature selection is added. 
The remaining features represent the input variables for the regression 
model. After the training, the selection mask for the features and the 
regression model parameters are fixed. Given a measurement from the 
test set, repeating the feature extraction, feature selection and regression 
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Figure 5.4 The principle of the FRM. 


steps with the fixed mask and the trained regression model parameters 
leads to the final estimate of the time-delay difference. 

A detailed description of the feature design, the feature selection ap- 
proaches, and the regression model follows in the next three sections. 


5.3 Feature design 


As the system identification in Chapter 3 has shown, the measurement sig- 
nals are band-limited and the contained information is time-dependent. 
Therefore, the signals can be sparsely represented in a different basis, e.g., 
the basis obtained by employing the AWPT. Since the AWPT coefficients 
represent a corresponding time-frequency window, they can be used 
to incorporate only relevant frequency parts and time ranges. In this 
section, three types of features are presented, two of which are used for 
the subsequent regression and one of which is an auxiliary variable used 
for the feature selection. Since the features are grouped in the form of 
matrices, the value range of the discrete frequency index k is changed 
from [0, K — 1] to [1, K] to be consistent with the matrix indices. 
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To begin with, the first step of the feature design is the AWPT of the 
reference and delay signals to get the corresponding coefficients 


c[k] = (s:n), ve lr]) > 
calk] = (salr] , vg nl), - (5.18b) 


(5.18a) 


Even though (5.18) implies that the coefficients are calculated via the 
inner product, the practical implementation is based on the dual-tree 
presented in Section 5.1.2, as the calculation using multirate filter banks 
reduces the computational effort significantly. To fully define the trans- 
form, both the filter types and the tree nodes must be specified. A prelim- 
inary study [151] has shown that, given the measurement signals with 
the excitation frequency 700 kHz and the sampling rate 50 MHz, the best 
results can be obtained by using nodes with a level s,, € [5, 6,7]. Fur- 
thermore, the study has shown good results if the filters of the first level 
(9o,ıp[r]) and the subsequent levels (g,,p[n]) are set to the Daubechies-20 
wavelets and the common factor solution of Selesnick with one vanishing 
moment and an approximation order of five, respectively. 

Using the calculated AWPT coefficients of the measurement signals, 
the features 


E(k] = |e, [k] — calf], (5.19a) 
=r c,[K] re 7 

Elk] = ee 1) Neel, (5.19b) 

Elk] = 5 (leskel + eal) (5.190) 


are calculated for each coefficient. Here, the coefficient index k repre- 
sents the time-frequency window—also known as the subspace—of the 
signal. In other words, after the signals have been transformed into the 
time-frequency representation, the three calculation rules (5.19) are ap- 
plied to each time-frequency window denoted by k. Note that taking the 
argument of the quotient of two complex values is equivalent to taking 
the difference of the arguments. The small constant e is added to the 
quotient to increase numerical stability on one hand and to improve the 
robustness against noise on the other hand, but it should be chosen well 
below the expected signal energy level of the target signals. Additionally, 
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note that the third feature (5.19c) is only used during the feature selec- 
tion and not for the regression, which is why in the following solely the 
features (5.19a), (5.19b) are examined in more detail. 

The motivation to use the mean absolute coefficient (5.19c) is that this 
feature represents the square root of the signal energy fraction within 
the respective time-frequency window, which can be used to decide 
whether the SNR is high enough to allow reliable information extraction. 
The other two types of features are chosen since they are approximately 
linearly correlated with the time-delay difference 7, which can be proven 
by the insertion of the narrowband assumption into the inner product. 
As (3.29) shows, the difference signal between the reference and delay 
signal is proportional to the time-delay difference 7 if r is sufficiently 
small. The projection of this difference signal onto the basis function, 


(saln] —seln], ven), 
(A,wor * COS (woltsn zu Ta)) 0 VE nl), 


A,WoT E (cos (woltsn == Ta)) ’ wr nl), 


X T, 


Q 


(5.20) 


Q 


is equivalent to the difference between the coefficients, since the transform 
is linear. Moreover, calculating the difference removes the symmetric 
interfering signals. The motivation to use the phase-difference feature 
is based upon the fact that a time delay of the target signal corresponds 
to a phase shift of the respective complex-valued coefficients. For this 
relation to hold, the basis function has to be analytic. A proof can be 
given by applying Parseval’s theorem to the inner product (5.18). If the 
target signal is a pure sinusoidal signal with circular frequency wọ, taking 
the time-delay difference leads to the phase-shifted coefficients 


(sin + for/2], YE [n]) ze sin], YE n)a. (6.21) 


It is easy to see that the two coefficients are out of phase with each 
other by exactly wor. The difference of the phases is consequently linear 
with respect to the time delay r. However, it must be noted that, firstly, 
the target signals are only approximately pure sinusoidal signals and, 
secondly, the interfering signals have not been taken into account here, 
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Figure 5.5 Dependency of the two proposed feature types with respect to the time delays 
Tm Of the different measurements. The curves are created by plotting the features (5.19a) 
and (5.19b) for the subspace with the highest energy over the time delays for data set 
Dim... Each point represents one measurement contained in the data set. 


which leads to a distortion of both the slope of the linear relationship 
and the linearity itself when working with real measurement data. 

A visual example of the relationship between the features and the time 
delay is shown in Figure 5.5. Here, the absolute- and phase-difference 
features of the subspace with the largest energy are shown for the simu- 
lated signals from the data set D,;,, 5. Due to the varying temperature, 
the absolute time delay 7, also varies, which influences the slope of the 
straight lines through the varying amplitude of the signals within the 
time window spanned by the corresponding basis function. Nevertheless, 
the graphs are in good agreement with the desired linear relationship 
between the time delay and the features. Furthermore, the linearity can 
be improved by restricting the temperature variation within the measure- 
ment package, which is used for training. It must be taken into account 
that the regression may then only provide good estimates in the process 
conditions range in which it was trained. A quantitative study on the 
generalizability will be given in Section 6.6.5. 

Finally, the absolute and phase-difference features for every subspace 
are concatenated into a single feature vector 


= Hi ER”? (5.22) 
to enable further analysis in a vector-matrix notation. 
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5.4 Feature selection 


According to Section 5.3, the number of features is two times the number 
of samples in the input signal since the basis transforms of the reference 
and delay signals yield K = N coefficients c, [k], cg [X], respectively. Using 
all features would not result in a stable regression that generalizes well. 
Furthermore, due to the sparsity of the signals in the AWPT basis, many 
coefficients are nearly zero, which in turn means that the corresponding 
features contain little information. Therefore, a feature selection has 
to be performed before the remaining features are further processed 
by the regression. This section describes in the first subsection several 
SotA methods for feature selection, which are implemented and will be 
tested in the algorithm for the TDE problem. In addition, a new iterative 
selection method is proposed—this method is characterized by the fact 
that multiple quality criteria for the features can be combined in an easy 
manner. 
Firstly, a few preparatory definitions are introduced: 


= Since one feature vector results from one of M measurements, the 
feature vectors from multiple measurements, indexed by m, can 
be grouped into the feature matrix 


Eee im E ROM. (5.23) 


= A feature selection is fully described by the set M containing the 
discrete indices k that have to be chosen from the feature vector. 
Consequently, the cardinality |M| determines the number of se- 
lected features K,. 


= In vector-matrix notation, the selection mask M € {0,1}%s*? 


yields the selected features from the feature matrix by the multi- 
plication 


an 
= M 
=, = 


a = (5.24) 
(K,,M) (K,2K) (2K,M) 


To this end, the selection mask is constructed from M by concate- 
nating K, = |M| rows, where each row consists of a unit vector 
eT € {0,1}'** with exactly one 1 at the k-th position. This also 
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ensures that there is no more than one 1 per column to avoid a mul- 
tiple selection of the same features. The resulting feature matrix 
=, contains the remaining features after the selection. 


= A selection is usually performed by specifying a threshold y or by 
taking the best K, features. For a more compact representation 
of the latter, Top x (9) C Q defines the operation, that selects the 
highest K, values from a set Q where each element N, € Q is a 
real number. 


= Since there are two types of features for the same subspace k, it 
can be distinguished whether to select each feature individually 
independent of the subspace or to specify subspace-specifically 
whether to select both or neither of the features of a subspace. 
Therefore, two selection constraints are considered in the following: 
Firstly, if the features are selected independently from the subspace, 
this is referred to as individual feature selection (IFS). Secondly, if 
only the subspaces and consequently both features of the respective 
subspaces are selected, this is referred to as subspace selection (SS). 
The number of selected subspaces is half the number of the selected 
features when using the SS approach. 


5.4.1 State-of-the-art feature selection methods 


Feature selection methods can be divided into three categories: filter 
methods, wrapper methods, and embedded methods [39]. In this thesis, 
representatives of all three types are presented and implemented to 
test their performance in feature-based regression of time delays. An 
overview of the evaluated methods with the corresponding references is 
given in Table 5.1. A more detailed description of the individual methods 
follows in the remainder of this subsection. 

Without modifications, the presented SotA feature selection methods 
perform the IFS. To adapt these methods to the SS, the quality measures 
of both features per subspace are averaged and a new threshold or top-N 
strategy is used to decide which subspaces to select. However, in addition 
to the constraints imposed on the selection, the averaging approach can 
lead to worse selection performances, due to the nonlinearity of some 
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quality measures. Moreover, some methods can only be used with the 
top-N strategy. 


Filter methods 


The filter methods are characterized by the fact that each feature is evalu- 
ated individually by a quality measure. A selection can then be obtained 
by setting a threshold y or by taking the features with the K, highest 
quality measures. In this context, the threshold or the number of features 
to select is a hyperparameter of the selection method. Since the subse- 
quently used regression model is not taken into account in the evaluation, 
the computational effort is generally lower than for wrapper and embed- 
ded methods. In this thesis, four filter methods are investigated for their 
applicability to feature-based regression of time-delay differences: the 
Pearson correlation coefficient (PCC), the F-test, the Laplacian score, and 
the RReliefF algorithm. Except for the Laplacian score, all these meth- 
ods are supervised, i.e., the labels 7,,, of the training set are used in the 
selection algorithm. 

The PCC—more precisely, the sample PCC—is the covariance between 
two variables normalized by their multiplied standard deviations. In 
other words, it is a measure of correlation between two quantities. There- 
fore, the PCC [98] 


= (Ekm = Em ~~ T) 
ET E Gn en 


Table 5.1 Overview of the used feature selection methods. 


Filter methods Wrapper methods Embedded methods 

= Pearson correlation » Wrapper with = Gaussian process 
coefficient (PCC) [98] linear regression regression (GPR) [96] 

a F-test [94] [62] a Least absolute shrink- 

= Laplacian Score [44] age and selection opera- 

= RReliefF [97] tor (LASSO) [119] 


= Neighborhood compo- 
nent analysis (NCA) [2] 
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is well suited to be a quality measure that ranks the importance of the 
features. In (5.25), the operator @,,, {-} denotes the sample standard de- 
viation with regard to m and Eu T denote the average feature and time 
delay over all M measurements, respectively. 

Another similar measure is given by the F-test, where a test statistic, 
which follows an F-distribution, is set up for each individual feature. In 
this selection method, a large F-value of the test statistic represents an 
important feature. According to Pirbazari et al. [94] the F-value for each 
feature, 


Jplk] = M2); (5.26) 


_ re lk] 


can be calculated based on its respective PCC r.,[k] and the number of 
samples M. If desired, the p-value can be derived from the cumulative 
distribution function of the F-distribution at Jp[k]. 

The only unsupervised method in the selection methods under inves- 
tigation, the Laplacian score, is based on a k-nearest neighbors (k-NN) 
similarity graph. Basically, for each feature vector, the k-NN feature vec- 
tors in | terms of the Euclidean distance are determined. Subsequently, 
& =E ill are inserted into a Gaussian 
kernel function to create the similarity matrix S = (s,,) € R“*™, with 

Sij = exp( ). The entries in the similarity ae that corre- 
spond to adonna ied B are set to zero. Finally, given the similarity 
matrix, the Laplacian scores are defined as 


Er > 1 Se 


M 7 
DI 1enS Im‘ Go 


Jlk]=1- (5.27) 


with the centered features € km: It has to be noted that a weighted arith- 
metic mean is used to center the features, where the weights are the row 
sums of the similarity matrix. A detailed description of Laplacian scores 
can be found in the paper by He etal. [44]. 

The last filter method, the RReliefF (regressional ReliefF) algorithm 
proposed by Robnik-Sikonja and Kononenko [97], is an extension of the 
original Relief algorithm to regression problems. It is based on rewarding 
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features whose values differ when the labels are different and punishing 
features whose values differ when the labels are similar. For each feature, 
the following steps are performed to get the quality measure. Firstly, a 
random observation is selected, and the k-NN in terms of the Euclidean 
distances of the feature vectors are determined. Secondly, the rewards and 
penalties of these k-NN are calculated and weighted by their respective 
Euclidean distance to update the quality measure. Lastly, these steps are 
repeated for further observations resulting in the final quality measure 
for each feature. More details about the algorithm can be found in the 
publication of Robnik-Sikonja and Kononenko [97]. 

By specifying the constraints of the selection—IFS or SS—and a thresh- 
old value y or the number of features to select K,, the selection mask M 
and its matrix counterpart M are obtained. In the case of SS, the selection 
matrix has a block diagonal matrix shape, i.e., the second half of the lower 
rows with the row indices K,/2 +1... K, is identical to the first half of 
the upper rows with the indices 1 ... K,. This is due to the concatenation 
of the features to a total feature vector. 


Wrapper methods 


Contrary to the filter methods, where the features are evaluated indi- 
vidually without considering the regression model, wrapper methods 
iteratively select subsets of the features and train the regression model. 
Consequently, the subset of features resulting in the least regression error 
is returned as optimal feature selection. However, since the amount of 
all subsets is the same as the cardinality of the power set, evaluating 
all subsets would require too much computational effort. Therefore, a 
greedy approach is usually used, where the features either are added or 
removed incrementally or a combination of both. In the literature, these 
approaches are called forward selection, backward elimination or step-wise 
selection [62]. Due to the greedy nature, the selection in each step is al- 
ways only locally optimal and the final selection does not represent the 
globally best solution. Furthermore, the chances of overfitting are high. 

As an example, the forward selection applied in this thesis is briefly de- 
scribed here. The algorithm starts with an empty set of selected features. 
Subsequently, for each feature, a simple regression model with only one 
independent variable is trained. The feature with the best regression 
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performance in terms of MSE or other performance criteria such as R? 
is added to the selected feature subset. These steps are repeated, where 
the trained regression model has an incrementally increasing number of 
independent variables, until a preset number of features is selected. Even 
if the number of regression models to be trained is greatly reduced in this 
way compared to the evaluation of all possible subsets, the computational 
effort is usually still significantly greater than with filter methods, espe- 
cially when a large number of features is to be selected. Due to the linear 
relationship of the features with the response variable (the time-delay 
difference), the linear regression model is chosen for the implementation 
of the sequential forward selection. 

Since the selection procedure is iterative and differs from the quality- 
measure-based filter methods, the SS cannot be obtained by averaging 
the quality measures. Therefore, the wrapper method is performed sepa- 
rately for both feature types and then, via a union and intersection of both 
selection sets, a subspace ranking is generated from which the subspaces 
are then selected via the top-N strategy. Threshold-based selection is not 
possible by means of wrapper methods. 


Embedded methods 


The last group of selection methods, the embedded methods, perform 
the regression training with a built-in feature selection method. From a 
computational point of view, they rank between the filter and wrapper 
methods. Again, the regression model is taken into account jointly with 
the feature selection, but it is also possible to specify a model for the 
feature selection and then discard it and use a different model for the 
actual regression. The three different examined embedded methods are 
the LASSO, the neighborhood component analysis (NCA), which was 
adapted to regressions, and the Gaussian process regression (GPR). 

Generally, the regression model parameters are penalized if they get 
too large. This restricts the regression model so that it is less prone to 
overfitting. A simple approach is the combination of a linear LS regression 
with a regularization leading to the objective of LASSO [119] 


M 


{ä,ö} = argmin (# > (m -0- Ta) + Az, D , (5.28) 


m=1 
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where the Lagrange multiplier \;, ensures that |lal|, < y for some thresh- 
old y. Higher Lagrange multipliers reduce the threshold resulting from 
the optimization, and can, therefore, be used to adjust the method to 
select fewer features, since more parameters a, are close to or directly at 
zero. An optimal Lagrange multiplier can be found by cross-validation, 
but to allow the number of selected features K, to be given as a hyper- 
parameter, the highest possible Lagrange multiplier is chosen such that 
the number of non-zero parameters is still greater than the specified 
number of features to select K,. Because the magnitude of the model 
parameters a, is dependent on the data set and the Lagrange multiplier, 
the LASSO-based feature selection can only be made via a top-N strategy. 
Furthermore, the SS-based feature selection is obtained by averaging the 
two normalized parameter sets that result when the LASSO is performed 
once each for the absolute and the phase features. 

The second embedded method is the NCA. Originally, it was proposed 
by Goldberger etal. [40] and applied as a feature selection method for 
classification by Yang etal. [135]. To make the NCA applicable to re- 
gression, the extension presented by Amankwaa-Kyeremeh etal. [2] is 
used in this thesis. The NCA is a non-parametric algorithm that aims 
at minimizing the average leave-one-out performance by learning the 
feature weights a,,, i.e., the average quality, if the label of a feature vector 
is predicted by consensus of its k NN, is maximized. Due to the length of 
the derivation and since the approach is not the focus of this work, only 
the objective function is stated here. The NCA finds the feature weights 
by the optimization [135] 


M M 
k 1 5 
=argmin— J Y (a) - (T, =T, À 3 
. nBr M m=1 j=l (Pm; (a) (Tm T3) ) T L lalı a 6 9) 


where p,,, ;(a) is the probability that the feature vector €; is chosen as 
the reference point for the label prediction of the feature vector €,,,. For 
m # j, it is calculated by 


X 2K 
exp (-0 i eet a? bem > Enzl) 


a == (5.30) 
yet jdm P (01 Dope alnm — Engl) 
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The weights a resulting from the optimization (5.29) represent the quality 
measures and can again be used to select the features by thresholding 
or by the top-N strategy. The SS-based selection is implemented in the 
same way as for the LASSO method. 

Finally, the last embedded method used is the GPR. In order to fit the 
Gaussian process to the training feature vectors and the corresponding 
labels, three practical assumptions are made: stationarity of the variance, 
the zero expectation value for all features, and a squared exponential 
covariance function with feature-specific variances as kernel [96]. The 
learned length scales of the automatic relevance detection kernel are 
used as quality measures of the respective features. 


5.4.2 Iterative subspace selection 


In this section, a novel iterative method for SS-based feature selection is 
presented, which is motivated by the fact that the information content in 
a subspace affects both feature types. The new method—in the following 
called the iterative subspace selection (ISS)—combines several problem- 
adapted quality measures in an iterative manner. The presented method 
uses only unsupervised quality measures and, thus, allows for a feature 
selection without ground-truth labels. 

Firstly, the quality measures are presented. Due to the design of linear 
features and the transform into a basis where the signals are sparse, three 
objectives stand out: a good SNR, as linear as possible features without 
distortions, and high sensitivity. Therefore, the three quality measures 


M (&glk,m]-£,[R])(Esik,m]-&£;[k]) 
msi FmilEglk,M]}-m {Eslk,m]} E 
M : 
Jilk] = $ En- alk m] i=2, (5.31) 
52, {&,[k,m]} i=3 


are introduced, with the features averaged over the measurements € glk], 
£,[k] and the sample variance 57, {-}. Note that even though the mean 
absolute coefficients £,[k, m] (see (5.19c)) are not used as features for 
regression, they can be well used here to select subspaces with high SNR. 
In order to keep the indices of the features short, the feature iteration 
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indices are denoted in brackets and not in matrix notation. The first qual- 
ity measure J, [k] is the PCC between the absolute and phase-difference 
features £; and £4, where the measurements denoted by the index m are 
considered the observations. Due to the missing labels in unsupervised 
approaches, the linearity between the features and the labels cannot be 
assessed directly by calculating the PCC between the features and the la- 
bels. However, the distortions that can occur in the features, if interfering 
signals are present, differ in nature for the absolute and phase-difference 
feature types, since the absolute value follows a root function and the 
phase follows an arctangent. Due to the suppression of interfering sig- 
nals in the absolute-difference features and their linear correlation with 
the labels, it is assumed that a strong linearity between €; and € , also 
indicates a good corresponding subspace. The evaluations of the ISS 
later in Section 6.6.2 will show the usability of this approach. The second 
quality measure J,|k] represents the average absolute coefficient values 
of a subspace over all measurement signals. Due to the orthogonality of 
the AWPT, the absolute value of a coefficient can be interpreted as the 
root of the signal energy in the corresponding subspace. Therefore, this 
quality measure is ideal to assess the SNR, since the noise is considered 
to be uniformly distributed over all subspaces due to the assumption 
as AWGN. Lastly, the third quality measure J3[k] is the variance of the 
phase-difference signal. This measure has been chosen because the phase 
of the coefficients is independent of the signal energy and a high value 
indicates a high sensitivity of the subspace to the measurement effect. 
However, the variance of the phase-difference features is very high in sub- 
spaces with low SNR. This and the strong sensitivity to signal distortions 
make this quality measure the least robust. 

Figure 5.6 shows an example of the three quality measures resulting 
from the measurement data set used in Section 3.3.1 for the identification 
of the measurement system. In order to visually improve the quality 
maps, they are displayed in logarithmic scale. Furthermore, the sub- 
spaces are ordered according to their time-frequency representation. 
From this example, several conclusions can be drawn, which emphasize 
the presumed properties. 

Firstly, the inter-feature PCC (Fig. 5.6 on the left) has a hard edge, be- 
ginning exactly at the absolute time delay r, at which the target signals 
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Figure 5.6 Normalized quality maps for time-frequency coefficients (in logarithmic scale). 
Left: PCC between phase-shift feature and absolute-difference feature. Middle: averaged 
absolute coefficient values. Right: variance of the phase-shift feature. 


of the direct propagation path are expected. In the time-frequency repre- 
sentation, another hard edge can be observed at about 12 MHz, which is 
caused by the constant e in the phase-difference calculation (5.19b). This 
can be explained by the fact, that above 12 MHz the coefficient values in 
the subspaces are lower than the constant, which means that no correla- 
tion can be measured anymore and the PCC tends to zero. However, apart 
from this, a distinction between the quality of the subspaces contained 
within this boundary is hardly possible. The remaining speckle noise 
is a display effect resulting from the use of the logarithmic scale. The 
absolute values (Fig. 5.6 in the middle) indicate multiple higher-order 
harmonics, resulting from the frequency responses of the basis functions. 
Additionally, at about 30 1s—the time delay of the interfering signals 
due to the Lamb waves—an increase of the coefficient energy can be 
observed. Lastly, the variance of the phase-difference feature (see Fig. 5.6 
on the right) shows good agreement with the variance becoming high 
exactly when a target signal is present, e.g., between 60 us and 80 ps, or 
when only noise is present, e.g., in the time range 0 — 30 ps. 

In order to combine the advantages of the different quality measures 
and solve the robustness problem of the variance measure õ?, 1&,[k,m]}, 
the ISS is proposed. Compared to a weighted mean, the iterative approach 
allows more flexibility such as combining the top-N strategy with a 
thresholding. Moreover, the approach is generally easier to parameterize. 
As the name implies, the quality measures are used iteratively to select a 
subset M“) of the K subspaces in each iteration step, which then gets 
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further reduced in the next iteration. In total, there are three iteration 
steps specified by the rule 


MO ={keEN|1<k< K}, (5.32a) 
M® = {kEeM® | k] >n}, (5.32b) 
M® = {k EM) | Jo[k] > 1}, (5.32c) 
M® = {k € M® | Js[k] € Topy ({Islk} remo) } - (5.32d) 


The first selection is based on the most robust but least discriminating 
quality measure, the inter-feature PCC. The necessary threshold param- 
eter 7, for this can be calculated by sorting the quality measures and 
taking the point of the strongest curvature—the same approach as for the 
automatic threshold estimation in Section 4.4. After this step, all signal 
subspaces containing only interfering signals are removed. The next se- 
lection step aims to separate the subspaces with low SNR. To this end, the 
VISU shrink threshold for denoising proposed by Donoho and Johnstone 
[27] is applied to identify all subspaces that are considered as noise. The 
AWGN standard deviation is estimated by the median absolute devi- 
ation of a measurement signal with only noise contained. In practice, 
this can be realized either by recording a signal without excitation or by 
observing subspaces known to contain only noise, such as the beginning 
of the measurement signals. Finally, the last step sorts the remaining 
subspaces and takes the top K, subspaces. Since subspaces with only 
noise or interfering signals are already removed, the robustness problem 
of this quality measure is mitigated. It is assumed that any remaining 
interfering signal components in the subspaces are reduced due to the 
feature design or can cancel out by using multiple subspaces. Neverthe- 
less, this is largely dependent on the number of features K, to select, 
which is at the same time the hyperparameter to choose most carefully. 

A visual example of the iterative selection process is given in Fig. 5.7. A 
reduction from the initial 7500 subspaces—determined by the 7500 dis- 
crete time samples—to the final 50 subspaces is shown from left to right. 
The negative PCC values of the first quality measure are a numeric effect 
and can be considered physically illogical. They are therefore removed by 
the threshold value y,, which has been found by the automatic algorithm 
here to be approximately 0.95. The final time-frequency mask shows 
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Figure 5.7 Iterative selection of coefficients. From left to right: linearity between features, 
absolute value of coefficients, variance of phase-shift feature. Top row: quality maps. Middle 
row: sorted and squeezed quality map. Bottom row: resulting mask after selection step. 


are good agreement with the measurement signal composition, since 
the commonly selected frequency is near to the excitation frequency 
(700 kHz) and three accumulation points can be observed in the time 
axis, which can be well explained by the three propagation paths of the 
ultrasonic signals. 

Since only subspaces are selected by the ISS, a post-processing is neces- 
sary to get the selection mask M. The final selection mask M™) resulting 
from the ISS may be denoted in matrix notation as follows: 


m 


M = [e,,...,€% nl" € {0,1} 2X. (5.33) 
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Because the ISS selects from subspaces and the total number of selected 
features K, should be consistent with the other selection approaches, 
the dimensionality is (K,/2, K). As already mentioned, the entire selec- 
tion mask M is obtained by taking both feature types for the selected 
subspaces resulting in 


Mo 
M = - | € {0,11% 2K, 5.34 
ROET en 
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The last component of the FRM is the regression model. As the features 
are designed to be linear with respect to the time-delay differences, linear 
regression models are chosen. Here, the features selected with the proce- 
dures described in the previous section represent the input variables. In 
addition, just as with the selection procedures, another alternative model, 
the GPR, is introduced. Subsequently, different training strategies to de- 
termine the parameters of the regression models are derived. Finally, an 
optimization-based approach is introduced to obtain the corresponding 
labels in order to use the FRM in an unsupervised fashion. 


5.5.1 Regression model 


In regression problems, a high-dimensional input vector, which is also 
known as independent variables or predictors, has to be mapped to the 
labels such that the estimation error is minimal. In the present problem, 
the input vector corresponds to the feature vector and the labels are the 
time-delay differences. Let the feature vector s m € R* be the vector 
that results from the feature selection of the m-th measurement. Then 
the three regression models 


Fm = Esma +o, (5.35a) 
Tm = Esma, (5.35b) 
Fim = GPR{Esmi Es TH (5.35c) 


allow an estimation of the time-delay difference 7,,, of the m-th measure- 
ment. While the former two models are both linear—the only difference 
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being the offset—the latter is a more complex model that can be best 
described as using the expectation value of a multivariate Gaussian 
distribution for the estimation. The Gaussian distribution is thereby de- 
termined by the given feature matrix =, and the labels 7 of the training 
data set. 

Physically, the model (5.35b) makes the most sense, since the features 
are designed to be zero for T = 0. However, this is only valid in the case 
without noise. Especially, the absolute-difference feature €; shows that 
even if the target signals and the interfering signals are suppressed due 
to the difference operation, the feature still contains the signal energies 
of the noise signals n, (t), nı(t) projected onto the respective subspace. 
Even though only subspaces with a high SNR should be selected, they 
still contain noise because the noise signals are distributed uniformly 
across all subspaces. Therefore, adding an offset may improve the over- 
all performance across the entire range of possible time-delay differ- 
ences. Nevertheless, the linear models are based on the assumption of 
linearly correlated features, which in turn depend on the assumption that 
T X 1/ fo and on the neglect of noise. In summary, the linear models may 
be too restrictive and are expected to perform better in the center of the 
time delay range covered by the training data than near the boundaries, 
but at the same time due to the restrictive model the risk of overfitting of 
the regression is reduced to a minimum. 

Since the linear models may be too restrictive, the GPR is also investi- 
gated and compared to the linear models. It has more degrees of freedom 
and also allows for an embedded feature importance ranking—called 
the automatic relevance detection—if the different features are weighted 
differently in the kernel function. However, training a GPR requires much 
more computational effort, which is why the model is combined with a 
preceding feature selection. 
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5.5.2 Training strategies 


The parameters of the used regression models have to be estimated 
during training. For the linear models, three parameterization methods 
are proposed, all based on a variant of 


T = er ’ a +0: lan +n. (5.36) 
(M,1) (MK) (Ku) (may (M,1) 


The first approach is to simply solve (5.36) globally via the LS estimator 


met 51 lo 
a =, stm =; 
= T >T T T. (5.37) 
fl Fee M | | 


It follows, of course, that this method of parameter estimation has the 
lowest mean squared deviations on the training data set. In the second ap- 
proach, the auxiliary variables @,, and 0, are introduced, which represent 
the slopes and the offsets if each feature is considered individually. Sub- 
sequently, the individual features € are mapped onto the time-delay 
differences domain by 


s,km 


T= (Ten B= meiste With Tim, = By ° Esem + Ök- (5.38) 


Finally, the individual slopes a, and offsets 0, are weighted to get the 
final model parameters 


a=woa, (5.39a) 
o=wTö, (5.39b) 


using the Hadamard product © and the weight vector w = v,/ ||v,|,, 
with v, € R” being the first eigenvector of the eigenvalue decompo- 
sition of the covariance matrix cov(T) = VAVT. This is motivated by 
the fact that the PCA finds the direction in which the most variance of 
the data can be observed. Therefore, less important features are automat- 
ically given a lower weighting factor. The last open question is how to 
find the individual slopes and offsets. One approach is once again the 
LS estimation, when only the observations of a single feature, denoted 
by the feature row vector & = R!*™ are considered. This leads to K f 
separate LS estimations 


a T —1 T 
ay = Es eE, Espl mi Sk oF, (5.40) 
Ok 11&s,k M Tin 
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Alternatively, the slopes and the offsets can be estimated using the 
median-based approach 


Teer 
a, = median,, | — —— >, (5.41a) 
Eskin a é [k] 


0, =T—4G,- Esk? (5.41b) 


which is more robust to outliers than the LS-based estimation. In (5.41), 
the average of the time-delay differences and the features over the mea- 
surements are denoted by 7 and E p, respectively. Both the median- and 
the LS-based estimation of the slopes a, can also be adapted to the 
more restricted model without offset (5.35b). To do this, simply remove 
from equation (5.40) the columns and rows containing the all-one matrix 
1,7 and from (5.41) the average values 7, E, „. In the following, these 
weighted regression training strategies are abbreviated as the weighted 
least squares (WLS) regression and the weighted median (WM) regres- 
sion. 

Compared to the linear models, the GPR model (5.35c) is harder to 
train. The parameters that have to be estimated are the kernel function 
and the noise variance. As fitting a Gaussian process to training data is 
not the focus of this thesis, only a short description of the assumptions 
and the optimization are given here. For amore detailed explanation, see 
the work of Rasmussen and Williams [96]. Firstly, no prior is given, i.e., 
the expectation values of the weights in the Bayesian linear regression 
model are set to zero and the variance is assumed to be stationary—this is 
also known as simple kriging. Secondly, the squared exponential kernel 
is used as the covariance kernel function. Based on these assumptions, 
the log marginal likelihood is minimized using a quasi-Newton method 
with the Symmetric Rank 1 method to approximate the Hessian matrix. 
Cross-validation is not used as the training takes several hours for only 
one run. 


5.5.3 Unsupervised time delay estimation 


In practical applications, the requirement of having a training data set 
with corresponding labels is a hard constraint that severely limits the 
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applicability of the previously discussed methods. However, there is no 
known method to train a regression model without access to the labels 
of the training data set. Therefore, another approach for the TDE, which 
is based on the optimization of an extended AMDF, is presented. The 
aim is to estimate the time-delay differences solely from the selected 
feature matrix =,, which is why, in principle, this approach could be 
used without the FRM. Nevertheless, the approach presented in the 
following is only used to generate the labels to the training data set, 
which are subsequently used for the training of the regression models. 
There are two main reasons for this. First, solving the optimization prob- 
lem package-wise limits the dynamics of the measurement system in 
an online—also called real-time—application and requires much more 
computational effort than the simpler linear regression. Second, while, 
through the regression-based method, the TDE can be performed on 
a single measurement consisting of a delay and reference signal, the 
optimization approach only allows the estimation of multiple time-delay 
differences for an entire measurement package. 

In the following, the objective function to minimize is derived. The 
idea is to get an estimate of the target signals from both the reference 
and the delay signals 


8¢(nt,) = salnt; +7,,/2;m) — §,(nt, +7,,/2), (5.42a) 


m 2 


Slnts) =s, (nt, = 727 out) — Sint TR), (5.42b) 
which are obtained by removing the estimated interfering signals 5,(nt,) 
followed by a time shift with the estimated time-delay difference 7,,,. 
As these two estimated target signals should be identical, if the interfer- 
ing signals and the time-delay differences are estimated correctly, the 


objective function to minimize can be set to 


N-1 
J (7,3;{nt,)) =D, [88 (nt,) — H(nt,)|, (5.43) 
n=0 


which can be interpreted as the time-discrete AMDF with the input sig- 
nals freed from interfering signals beforehand (cf. (2.9)). However, using 
the objective (5.43) has several drawbacks. As the signal model (3.25) 
has already shown, the time-delay difference differs for different time 


117 


5 Feature-Based Regression Method 


ranges, i.e., the model is not globally valid. Furthermore, the frequency 
characteristics of the target signals are not taken into account and the 
optimization is too high-dimensional to converge due to each interfering 
signals sample 8,[n] being a parameter to optimize. These drawbacks are 
alleviated by the adaption that the objective is solved in the subspace 
spanned by the best basis function wf [n]. Thus the objective function is 
reformulated as 


= |) 5d{n]C[n] — s: [nC In] 


n 


=0 
(5l), wen), - (el, ven), |. 64 


Similar to before, the inner products of the basis function with the input 
signals represent the AWPT coefficients c4[k], c,[k] of the respective 
subspace. Since the basis function has compact support in the time- 
frequency range, only the frequency and time range of the target signals 
are considered by the objective function (5.44). Moreover, due to the 
orthogonality of the AWPT, only interfering signals of the form 


5,[n] = gyp. [n] + mym in] (5.45) 


result in non-zero coefficients if the problem is projected onto the sub- 
space. This reduces the degrees of freedom of the interfering signals to 
two. In order to estimate the time shift 7,, with subsample precision, 
the inner products are evaluated via Parseval’s theorem in the Fourier 
domain. Inserting (5.42) and (5.45) into the objective function (5.44) and 
transforming the result into the Fourier domain leads to 


Nk a 


Fn a =| (Sam), Weine") 
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with the time-shift-dependent factors 


n 
Ree) = (Vf ‚—2j We [ny] sin(n 2) , (5.47a) 


KG.) = (Wirt) -2 WE sin(n a) . (5.47b) 


In these equations, capital letters represent the DFT of the respective 
signals. A closer look at the objective function shows that three unknown 
variables have to be solved while only one measurement is available, 
rendering the problem underdetermined. 

Due to this, the objective function for multiple measurements are 
combined. The stationary-interfering-signals-assumption means that the 
factors er eo are constant for different measurements—assumed that 
the process conditions do not vary to a large extent. However, while the 
interfering signals are constant, the time-delay difference differs from 
measurement to measurement. A solution to this problem can be ob- 
tained from the PCA of the part of the feature matrix =, that contains only 
the absolute-difference features £5, since these features are unaffected 
by interfering signals. Even though their sensitivity is dependent on the 
signal amplitude in the time range of the corresponding basis function, 
the first eigenvector of the PCA v; = (Vm)m=1...m İS a good approxima- 
tion of the relative relationships between the time-delay differences of 
different measurements. 

Finally, the derivation of the objective function leads to the three- 
dimensional minimization problem 


A Ivila t 
TCE m= 1112 


M Tom 
Giese") = argmin XJ Fr hia) f (5.48) 
1 


which can be solved by the interior point algorithm. Consequently, the 
labels of the entire training data set can be obtained by 7,,, = Tp) Um / ||Villo- 
These labels are then used to train the regression models with the pre- 
viously presented methods. If the chosen selection method works in an 
unsupervised fashion and the labels are estimated through the optimiza- 
tion approach, the entire FRM is also unsupervised, rendering this class 
of TDE methods directly comparable to the SDMs and the SotA TDE 
methods. 
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The methods presented in the chapters 4 and 5 are based exclusively 
on the signal model and can be applied to a wide range of ultrasonic 
systems, where the objective is to get a high-precision estimate of very 
small time delays and the interfering signals are the main influence that 
limit the precision. The only constraint is that the assumptions on the 
stationarity of the interfering signals must hold and the target signals 
must show a varying time delay, when an entire measurement package 
is observed. 

In order to evaluate the proposed methods on experimental data, the 
following presents how the methods perform in an UFM setup. To this 
end, the measurement setup with six different data sets and the eval- 
uation metrics are described in the first three sections. Subsequently, 
the performances of the SotA methods are evaluated and the methods 
that perform best on the respective evaluation metrics will be used as 
baseline for the new methods—SDM and FRM. After the novel methods 
are evaluated in terms of systematic and stochastic errors for different 
hyper parameter constellations, the chapter is concluded with an overall 
discussion and comparison of the different methods. 


6.1 Measurement setup 


For the UFM system, eight ultrasonic clamp-on transducers (UT1,...,UT8) 
are connected to a pipe, which was made of 4mm stainless steel with 
8cm diameter. The transducers are arranged in two circles, i.e., four 
transducers per circle, with the circles being 6 cm apart in axial direction. 
Four transducer pairs (TP1,...,TP4), which show a measuring effect, i.e., 
Arx # 0, can be built with this configuration. In order to obtain a data set 
with a physical damping of the interfering signals, another identical pipe 
with the same transducer configuration but an additional damping mat 
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Figure 6.1 Setup of the UFM system. For better illustration, in axial direction (left) the 
cross-sectional view is given and in radial direction (middle) only the transducers that 
build a TP are drawn, i.e., the ultrasonic transducers 3, 4, 7, 8 are not visible in the axial 
cross section and the ultrasonic transducers 5, 6, 7, 8 are not visible in the radial view. 


mm 
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is constructed. As damping material aluminum butyl—a self-adhesive 
butyl rubber with an aluminum carrier foil—is used. The transducer and 
pipe setup with an overview of the dimensions is depicted in Fig. 6.1. For 
better illustration, only four of the eight transducers are shown in both 
the axial and radial view. Additionally, possible ultrasonic propagation 
paths are included in the setup, e.g., the direct path (violet), an axial 
reflection (dashed green), and a radial reflection (green). Note that the 
interfering signals propagate through the pipe wall (see the dashed red 
path in the radial view) and that the direct path’s radiation angle into 
the medium is subjected to Snell’s law. 

The pipe prepared in this way, also called measurement section, was 
inserted into a circular flow with a pump to control the flow, a heating 
element to control the temperature and a reference flow meter to get the 
ground truth for the VoF. By installing a flow straightener in combina- 
tion with a sufficient inlet length right before the measurement section, 
a fully developed turbulent flow profile was achieved. This system is 
filled with water or ethanol, depending on the data set to be recorded, 
and any remaining air is removed as best as possible via a bubble sepa- 
rator. Sometimes, especially at the beginning of the experiments, a few 
remaining stuck bubbles may come loose during the experiment, which 
distorts the measurement signals due to the scattering of the ultrasonic 
waves. The measurements recorded during these incidents are removed 
before the evaluation, which is why sometimes a few data samples are 
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Figure 6.2 Typical measurement signal with highlighted signal composition. 


missing. In terms of temperature, the experiments are set in such a way 
that the measurements go up to about 30 °C, since the overall system can 
then be considered safe, even if, for example, ethanol is used as medium. 
As starting temperature, the lowest possible value is chosen, which is 
dependent on the room temperature since no cooling element was in- 
stalled. The exact temperature curve during the experiments depends 
on the specific heat capacity of the medium, since the heating element 
was not installed in closed-loop control. 

The data acquisition was performed at a sampling rate of f, = 50 MHz 
using a preamplifier and a PXIe-1062 station with a PXIe-5171 ADC 
module and a PXI-5412 DAC module, which also controlled the pump 
and the heating element. Combining all controls, the excitation, and 
the recording of the ultrasonic waves into a single system allowed hav- 
ing synchronized measurement readings and reduced jitter. In order 
to get the ultrasonic recordings of the downstream (reference signal) 
and upstream (delay signal) direction as simultaneously as possible, 
and in addition, to perform measurements with multiple TPs under the 
same process conditions, the excitation and recording order is as fol- 
lows: UT1/UT2; UT2/UT1; UT3/UT4; UT4/UT3; UT5/UT6; UT6/UT5; 
UT7/UT8; UT8/UT7. In this notation, the first transducer is operated 
as sender and the second transducer as receiver, respectively. Since the 
reverberation of the ultrasound in the pipe should not affect the next 
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measurement, the time between two consecutive measurements is set to 
10 ms even though the absolute recorded time of a single measurement is 
only 150 ps. An example measurement is depicted in Fig. 6.2, where the 
estimated composition of the measured signal is highlighted. Note that 
this highlighted composition is only manually estimated based on the 
signal envelope, the geometric relationships, and the present time-delay 
differences. As for all ultrasonic measurements in this thesis, the sending 
transducer was electrically excited using (3.18), with the center frequency 
being 700 kHz and the bandwidth being 245 kHz. 


6.2 Data sets 


Six different data sets were recorded to cover a wide range of scenarios. 
The data sets can be classified according to whether the pump was set to 
create a constant flow or a varying flow. Furthermore, data sets with two 
different media—water and ethanol—are recorded, once in an undamped 
pipe and once in a damped pipe. An overview of the data sets can be 
found in Table 6.1. The separation of the data sets in constant and varying 
VoF is important since the FRM requires a variable time-delay difference 
in the training data set, which is why only data sets with varying VoF are 
later used to train the FRM. The media influence on one hand the absolute 
time delay 7, due to their different SoSs and on the other hand the 
acoustic load of the ultrasonic system, causing the level of the interfering 
signals to change. The last distinguishing point of the data sets, the 
damping mat, reduces the interfering signals by a certain factor but also 


Table 6.1 Overview of the used data sets and their properties. 


Data set ID Medium Pipe Type of VoF 
Die water undamped v constant 
Di, water undamped v variations 
Do water damped v constant 
Do, water damped v variations 
D3 ethanol undamped v constant 
Day ethanol undamped v variations 
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Figure 6.3 Time delay and temperature process conditions of the recorded data sets. 


causes the remaining interfering signals to be no longer stationary since 
the attenuation of the selected material is temperature-dependent. 
Figure 6.3 shows the time-delay differences and the temperature values 
during the measurements. For the constant flow data sets, the VoF was 
set to 1ms ‘and for the varying flow data sets, the three flow velocity 
levels 0.3ms~', 0.6ms ', 1.2ms”! were repeated cyclically. Since the 
experiments were terminated at about 30°C and the different media 
and pump control strategies resulted in temperature rises of different 
degrees, the number of measurements in the data sets is not consistent. 
Furthermore, the time-delay difference varies even in the case of constant 
flow, because varying temperature leads to varying SoS and this, in turn, 
affects the time-delay difference, cf. Eq. (3.26b). Especially for the medium 
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ethanol (D3,, Dy), the temperature increases rapidly due to the lower 
specific heat capacity, which is why the data sets for ethanol contain only 
less than 1500 measurements, while the data sets for water contain more 
than 2000 measurements. Also the time-delay differences of the ethanol 
data set are generally larger due to the smaller SoS. 

Contrary to the ethanol data sets D,. and D3,, the water data sets with 
the damping mat attached to the pipe Dz, and Da, do not differ in terms 
of process conditions from the water data sets without damping mat. 
The only difference is the level of interfering signals and the fact that the 
assumption on their stationarity is not met. 

The last point to mention is that for each measurement in the data sets 
there are ultrasonic measurements from four independent TPs, each of 
which can be considered as a realization of the stochastic process that 
determines the shape of the ultrasonic signals. Here, the exact shape of 
the ultrasonic signals is considered random because they are sensitive 
to minute variations in material parameters, acoustic coupling, and the 
geometric arrangement of the transducers. Therefore, the probability to 
obtain exactly the same ultrasonic signals for two different TPs is nearly 
zero, even if the same transducers are removed from the pipe once and 
reattached to the pipe. However, the basic properties of the ultrasonic 
signals, such as the time of arrival of the target signals or the time-delay 
difference, are preserved. 

Finally, the notation of the data sets is explained: each used data set 
is completely defined by the notation D,.ırpı, where the symbol before 
the vertical bar denotes which data set and the number after the vertical 
bar denotes which TP to use. In this example, the first TP of the constant 
flow water data set in the undamped pipe is meant. 


6.3 Preliminaries and evaluation metrics 


The performances of the proposed TDE methods are evaluated in terms 
of accuracy, robustness against noise, and possible requirements for 
training data sets or necessary signal dynamics. For this purpose, several 
metrics are introduced to describe the performance. Furthermore, the 
performances are investigated while varying the available signal dynam- 
ics, training data sets, and SNR. Before the methods can be compared, 
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a few preliminary tasks have to be solved. In addition to determining 
the ground truth, this also includes estimating the time ranges to be 
evaluated for the time-based methods—which essentially include all 
SotA methods and the SDMs—or characteristic values such as the SNR 
or the CRLB. 

To begin with, the ground truth is to be determined. To this end, a 
reference measurement system from Endress+Hauser yields the aver- 
age VoF v, which has to be converted into a time-delay difference via 
reformulation of (3.4) to the time delay 7. Additionally, the hydraulic 
correction factor (3.8) has to be multiplied by the average VoF to get the 
average velocity of the flow along the ultrasonic transmission path. In 
summary, the ground truth is calculated by 


2k (v, T)-v- L- cos(a) 
T = 
c? (T) 


; (6.1) 


with the hydraulic correction factor [48] 
k (T, T) = 1.125 — 0.0115 - logy, (Re(0,T)) , (6.2) 
which in turn depends on the Reynolds number 


le) 
n(T) 


The geometric sizes, such as path length L and beam angle a in (6.1) 
are calculated by Snell’s law, the transducer placement, and the pipe 
geometry. The other fluid properties —the SoS, the density p(T) and the 
dynamic viscosity n(T)—are determined according to the NIST webbook 
[67]. Note that since the Reynolds number is Re > 10* for all temperatures 
and all VoFs during the experiments, the flow is always fully turbulent. 

Next, the expected effects of the interfering signals on the systematic 
error of the TDE and a good metric that reflects the interfering signal level, 
in this thesis called ScNR, are discussed. Once again, the narrowband 
assumption is used to derive the expected estimation errors. As already 
shown in the introduction, the additive interfering signals lead to a phase- 
dependent error of the TDE, if the narrowband assumption holds (cf. 
Fig. 1.4). However, these results are only valid for large variations of the 


Re(v,T) = (6.3) 
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Figure 6.4 Relative estimation error depending on the absolute time-delay difference for 
different ScNRs and time-delay differences. 


time delay, such as for the absolute time delay in the UFM application. 
When small time delays are present, such as the time-delay differences 
in UFM, the effects on the systematic error of the time-delay difference 
estimate can be described by calculating (1.4) once for the reference and 
the delay signal and then taking the difference between the results. This 
is equivalent to calculating the gradient of the curves shown in Fig. 1.4 
by the secant approximation. Mathematically, the estimation error can be 
expressed via the difference of the arctangent functions and simplified by 
a Taylor series expansion, which is terminated after the linear member, 
leading to 


Ca) or, (64) 


with a phase and ScNR-dependent constant C and the respective esti- 
mated absolute time delays 


e aa O 


+ cos(wolTa +7)) + A; cos(y) 


a 


Figure 6.4 shows a numerical evaluation of the estimation error derived 
in (6.4) for two different ScNRs, each with three different time-delay 
differences. There are three insights to take from this figure. Firstly, the 
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different time-delay differences lead to the same relative error, which jus- 
tifies the use of the percentage error as the evaluation measure. Secondly, 
the percentage error is dependent on the absolute time delay, which 
is why the evaluation of a TDE method in UFM may only be done by 
looking at the percentage errors over all absolute time delays covering 
one period 1/ fọ. Thirdly, even for very small interfering signal levels 
(ScNR = 40 dB), the estimation error amounts up to 1%. In summary, 
when considering the estimation error plot as a function of temperature, 
an oscillating error is expected, whose amplitude depends on the ScNR, 
since a varying temperature leads to a varying absolute time delay. 
Due to the estimation error being approximately proportional to the 
ground truth of the time-delay difference, the performance of the TDE 
methods is evaluated in terms of percentage error for each measurement 
Elm] = 100% - La (6.6) 
Tm 
Following the argumentation above, to get the performance of a method 
independently from the current phase between target and interfering 
signals, the mean absolute percentage error 


MAPE=— © |Eim]| (6.7) 
Ta, met, 
is calculated, with the absolute time delays of the test measurements 
being drawn uniformly from the interval 


to = [72,0 ’ T2,0 at 1/ fo) ’ (6.8) 


which contains approximately one period of the signals. In equations 
(6.7)-(6.8), the measurement with the smallest available absolute time 
delay in the data set is considered the beginning of the interval T, ọ and 
the number of measurements taken into account for the MAPE is M. A 
special case, that will be later used, is the evaluation of the FRM. Since 
it is a learning-based method, a distinction should be made between 
a global quality and a local quality when examining the performance. 
For the local quality, the MAPE includes only test measurements, whose 
absolute time delay is contained in the interval 

t= [min(r Tein), max(r} an] ; (6.9) 


a,m m 


Ta ‚m 
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which is spanned by the absolute time delays contained in the used 
training data. 

Additionally to the MAPE, which evaluates the estimation error on the 
entire test data set regardless of its origin, the standard deviation of the 
stochastic percentage error (STDSPE) is introduced as a metric to assess 
the stochastic error decoupled from the systematic error. To this end, after 
the systematic error is removed via subtraction of the lowpass filtered 
component, the STDSPE is defined as the sample standard deviation 


STDSPE = 5 „{Elm] — LP{E|m]}} . (6.10) 


Since the low-pass filtering should not result in a delay of the estimation 
error curve and should be robust against outliers, a moving median 
filter of order 41—where the odd order is advantageous to obtain a 
symmetric filter—was chosen. If this stochastic error is not expressed as 
a percentage, a theoretical lower bound, called the Cramér-Rao lower 
bound (CRLB), can be calculated. Although Carter [15] presented the 
equation to calculate the CRLB for TDE problems, the PSDs of the noise 
and target signals must be estimated with sufficient accuracy, which is 
not possible due to the measurement signals being too short. Therefore, 
the simplified formula according to Walker and Trahey [124] 


1+25NR 3 — 1 
SNR?  2mtoas 120, f2 +o? 4 


CRLB = (6.11) 


which assumes flat spectra with constant SNR and limited bandwidth, 
is applied to estimate the CRLB. Equation (6.11) is dependent on the 
observation time t „ps, the bandwidth o;, the center frequency f, and the 
SNR. Unfortunately, these quantities are also not directly known and 
must be estimated from the measurement signals. Consequently, the 
PSDs are approximated via Welch’s method using the same hyperparam- 
eters as for the investigation of the target signals’ spectral characteristics 
in Chapter 3. Based on the PSDs, the 3 dB bandwidth is obtained, i.e., 
the frequency range in which the PSD is above half its maximum value. 
Furthermore, the center frequency is chosen as the midpoint of the 3 dB 
frequency band. Since the AWGN is distributed over the entire time 
and frequency range, but the target signals are only present in a small 
time range and frequency band, the SNR is calculated as the ratio of 
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the signal energies in the considered time range and frequency band of 
the target signals. This has the advantage that only the part of the noise 
that interferes with the target signals is considered since the remaining 
part can theoretically be removed by filtering. It should be noted that 
the different methods each use a slightly different observation window 
and therefore the CRLB must be calculated separately for each method 
to obtain comparable results. 

The last preliminary task is to find the time range of the target signals. 
This time range is only necessary for the SotA methods and the SDMs, 
because the training-based FRM implicitly finds its useful time ranges 
through the feature selection. For simplicity and since the evaluation is 
based on the time-delay differences, whose ground truth is given by (6.1) 
only for the diametral path, only the direct propagation will be used. 
Because the algorithm to detect the time range needs to be very robust, 
several conditions are imposed on the algorithm to ensure that only 
meaningful time ranges are returned. However, a complete description 
of the algorithm would distract from the focus of the work, which is why 
only the basic approach is discussed here. Further details can be found 
in [157]. The basic principles of the algorithm can be summarized by the 
following steps: 


= First, the difference signal As[n] and the average signal of the 
reference and delay signals are formed. Basing the time range 
detection on the difference signal would be advantageous since the 
symmetrical interfering signals are eliminated, but the difference 
signal can only be used, if its signal energy is sufficient, i.e., the 
time-delay difference should not be too small. Otherwise, the time 
range detection is based on the average signal—the decision can 
be made through a threshold for the signal energy. Subsequently 
to the decision which signal to use, the envelope for that signal is 
calculated via a spline interpolation of the maxima, as this is more 
robust than applying the Hilbert transform. 


In the next step, all local maxima and minima of the envelope 
are determined. When detecting the extrema, threshold values 
must be set for the absolute height, the peak prominence, and the 
distance between two peaks in order to be robust against signal 
interferences or distortions. 
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a Finally, the local maximum, which is closest to the theoretical 
time delay calculated by the path length and the SoS, is taken 
to be the wave packet of the direct propagation path. Based on 
this rough determination, the upper limit of the time range is 
the next consecutive local minimum of the envelope, as this is 
considered the time, when the next wave packet starts. Either the 
local minimum before the determined local maximum, the time at 
which the envelope falls below a certain threshold or the minimum 
limit, which depends on the upper limit and the determined local 
maximum, is selected as the lower limit of the time range. The 
decision which lower limit to choose is made in such a way that 
the resulting time range is as small as possible. 


6.4 State-of-the-art methods 


In this section the baseline for the proposed TDE methods is determined. 
For each data set, the time-delay differences of the direct propagation 
path are estimated by using the time-based zero-crossing method (ZCM), 
and four different correlation-like methods, namely the cross-correlation 
(CC), the generalized cross-correlation (GCC) with the phase transform 
processor, the AMDF, and the ASDF. Since the time-delay difference has 
to be estimated with subsample precision and depends on the observed 
time range, the SotA methods are combined with the time range detec- 
tion presented in the previous section and with an interpolation. For 
the interpolation of the zero crossing, the method proposed by Kupnik 
et al. [63] is applied. In contrast, the CC, the GCC, and the ASDF use a 
quadratic interpolation around the maximum of the correlation function. 
Among the correlation-like methods, only the AMDF method is an excep- 
tion, since here the objective is to find a minimum of the AMDF, whose 
shape looks like the absolute value of a sine wave around its minimum. 
Therefore, a double linear interpolation is applied for the AMDF-based 
subsample TDE. 

The results of the ZCM, the CC and the AMDF for the four transducer 
pairs of data set D,. are shown in Fig.6.5. As already expected, the 
estimation error shows an oscillation with a period in the range 1— 1.5 ps, 
which shows a good agreement with the theory as a period duration of 
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Figure 6.5 Error curves of some selected SotA methods. The respective MAPE is given in 
brackets. Although the estimation error of all measurement signals contained in the data 
set D,, is shown here, the evaluation metric MAPE was only calculated on the set covering 
one period cycle according to (6.8). 


1/fo ~ 1.4 ps was predicted theoretically. The different period durations 
for the different transducer pairs can be explained by the slightly varying 
material characteristics of each transducer. Only the local frequency 
of the target signals and the interfering signals in the evaluated time 
range are relevant. Therefore, the local frequencies may deviate from the 
excitation frequency. Furthermore, this figure shows that even the level 
of the interfering signals may differ from transducer pair to transducer 
pair, since, e.g., the amplitude of the oscillating estimation error E[m] 
for TP4 is significantly smaller than for TP3. Apart from that, all SotA 
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Table 6.2 MAPE results of the different SotA TDE methods for all data sets. 


MAPE/% ZCM CC GCC AMDF ASDF 


TP1 6.1 5.6 5.3 52 5.6 

5 TP2 4.5 2.9 2.9 3.3 2.9 
u TP3 6.2 7.0 5.4 6.8 7.0 
TP4 2.0 2.6 17 2.5 2.6 
Average 4.7 4.5 3.8 4.5 4.5 

TP1 1.6 i4 1.2 1.3 1.1 

D TP2 1.3 1.2 0.8 1.9 1.2 
a TP3 1.8 2.1 1:3 2.3 2.1 
TP4 0.9 1.0 1.0 1.5 1.0 
Average 1.4 1.3 1.1 1.8 1.3 

TP1 183 18.7 183 18.7 18.7 

D TP2 7.7 6.9 5.5 8.1 6.9 
3c TP3 7.8 7.8 6.5 9.0 7.8 
TP4 13.0 13.9 12.7 13.9 13.9 
Average | 11.7 11.8 10.8 12.4 11.8 


methods investigated here perform essentially the same with only slight 
differences. While the ZCM is best in terms of MAPE for two of the 
transducer pairs, the AMDF and the CC are better for the other two 
transducer pairs. Normally, the zero-crossing-based TDE would be more 
noisy, but using the approach of Kupnik etal. [63], the information of 
a whole period can be taken into account to get a zero crossing that is 
more robust against noise. For this reason, no obvious difference in the 
stochastic estimation error can be observed. 

The MAPE results of the remaining SotA methods on all data sets are 
listed in Table 6.2. While the MAPEs of the respective best methods are 
highlighted in green, the average values over all transducer pairs are 
printed in bold. The best values of the respective data sets represent the 
baseline when, in the further course, the results of the SDMs or FRMs 
are given for individual transducer pairs or are averaged over an entire 
data set. It is noticeable that the damping mat has significantly reduced 
the interfering signals, because the MAPE values of the data set D,, are 
significantly lower compared to the data set without damping mat D,,. 
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Table 6.3 STDSPE results of the different SotA TDE methods for all data sets. 


STDSPE/% ZCM CC GCC AMDF ASDF 


TP1 0.60 0.58 0.59 0.56 0.58 

TP2 0.56 0.56 0.57 0.53 0.56 

7 TP3 0.59 0.57 0.59 0.51 0.57 
TP4 0.62 0.60 0.61 0.56 0.60 
Average | 0.59 0.58 0.59 0.54 0.58 


TP1 0.58 0.57 0.57 0.55 0.57 

TP2 0.57 0.56 0.56 0.53 0.56 

á TP3 0.59 0.58 0.58 0.55 0.58 
TP4 0.61 0.59 0.60 0.56 0.59 
Average | 0.59 0.58 0.58 0.55 0.58 


TP1 0.70 0.67 0.77 0.64 0.67 

TP2 0.89 0.83 0.86 0.83 0.83 

7 TP3 0.71 0.69 0.81 0.71 0.69 
TP4 138 1.55 1.41 1.22 1.55 
Average | 0.92 0.94 0.97 0.85 0.94 


Furthermore, the ethanol data set D3 , seems to have the highest level of 
interfering signals, which can be explained by the worse acoustic cou- 
pling between the pipe wall and the liquid. Once again large differences 
can be observed between the different transducer pairs of the same data 
set. In summary, all SotA methods are in the same order of magnitude, 
with the GCC method showing slight advantages on average for all data 
set. 

Complementary to the MAPE, which mainly evaluates the system- 
atic error—except in the case when the systematic error is significantly 
smaller than the stochastic error—the STDSPEs for all data sets as a 
quality measure for the stochastic errors are given in Table 6.3. Here, it 
can be seen that, as already suspected from the error curves in Fig. 6.5, 
the stochastic errors of the methods are very similar, with the AMDF 
being slightly better than the other methods. Interestingly, the difference- 
based correlation of the AMDF can handle noise better, but against the 
interfering signals, it performs rather averagely. 
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6.5 Signal-dynamics method 


The SDMs can be separated into the PCA-based, the B-spline(BS)-based 
and the joint-B-spline(JBS)-based point cloud processing, each of which 
can be further separated into the direct TDE estimation approach, the 
local interfering signals compensation and the global interfering signals 
compensation. The difference between the latter two is that in the lo- 
cal case, the interfering signals are compensated for each measurement 
package individually, while in the global case, all estimated interfering 
signals are combined, freed from outliers, and then averaged to get an 
estimate that is used for the interfering signals compensation of all mea- 
surement packages equally. Firstly, the direct TDE methods on a single 
data set are investigated, examining both the expected properties and the 
dependency on the process conditions. In the next step, the concept of 
compensating the interfering signals locally and globally is proven. Fur- 
thermore, visual examples of the estimated interfering signals, the outlier 
detection, and the globally estimated interfering signals are shown. This 
section concludes by examining how to specify the hyperparameters. To 
this end, the approaches are applied to the water data sets using all hy- 
perparameter combinations. Then, using the optimized hyperparameter 
set, the performances of the SDMs are averaged over all data sets and TPs 
and compared with the GCC method—the best SotA method in terms of 
MAPE. 

As in the rest of the thesis, for the SDMs to be applicable, the mea- 
surement signals of the data sets were divided into packages with a 
signal dynamic of 50 ns, i.e., within the package there are measurement 
signals whose absolute time delays cover an interval of the width 50 ns. 
An exception to this is the investigation of the influence of the signal 
dynamics in Section 6.5.3. Furthermore, the packages were selected from 
the entire data set with a variable overlap such that they uniformly cover 
the entire data set. 


6.5.1 Direct time delay estimation 


In the first evaluation, the direct TDE methods are applied to the water 
data set with constant flow D, .ırpı to check whether the approach 
can also be applied to real data. After the signals were separated into 


136 


6.5 Signal-dynamics method 


PCA-based BS-based JBS-based 
80 
| | | | Ground truth 

70 |- A F | ||—SDM 
a 
= 60 j} 4 F 4 
cb Nel 

50 f* | y | F | 

40 | | | | | | 

68.5 69 69.5 7068.5 69 69.5 70 685 69 69.5 70 


T,/ps T,/ps Ta/ns 


Figure 6.6 Results of the direct TDE using the SDM for the data set D,.rpı- The signal 
dynamics of the packages are set to Ar, = 50 ps. The median absolute time delay r, m of 
each signal package was distributed uniformly on the whole process condition range. 


packages and the relevant time ranges were extracted package-wise by 
the method described at the end of Sec. 6.3, the time-delay differences are 
estimated using the PCA-based, the BS-based, and the JBS-based point 
cloud processing, respectively. While for all methods the same signal 
dynamics (50 ns) were specified, the overlap for the PCA-based TDE was 
increased to get a densely sampled error curve comparable to the BS- and 
JBS-based SDMs. This is necessary since the PCA-based direct TDE can 
only estimate one average time-delay difference per package. The last 
hyperparameter to specify is the point cloud generation for the BS- and 
the JBS-based approaches. According to Equation (4.47), there are two 
possibilities: either the point clouds of the delay and reference signals are 
generated and processed separately, or the delay and reference signals 
are concatenated to get a combined point cloud. In this evaluation, the 
point clouds were concatenated. 

The results are depicted in Fig. 6.6. Although all three methods basi- 
cally estimate the time-delay difference well, there are slight differences, 
which can be explained by the different properties. The PCA-based SDM 
is less noisy. This is to be expected because the time-delay differences 
are averaged over the whole signal package, which can be interpreted as 
an individual estimation with subsequent moving average filtering. Fur- 
thermore, the joint estimation of the BSs for two consecutive time steps 
obviously also reduces the stochastic component of the TDE. However, 
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Figure 6.7 Results of the direct TDE using the SDM for the data set Dj) pp, with varying 
time delay. The signal dynamics of the packages are set to Ar, = 50 ps. The performance 
of the PCA-based SDM deteriorates significantly because the requirement—a constant 
time-delay difference—is not met at all. 


in terms of systematic error, the BS-based SDM shows better results than 
the PCA-based and the JBS-based SDMs. The larger systematic error of 
the PCA-based SDM can be explained by the fact that the constant time 
delay assumption does not perfectly hold, since the time-delay difference 
varies even in constant flow scenarios due to the varying SoS (see Fig. 6.3). 
Contrary to this, the JBS-based approach relies more on the assumption 
that the interfering signals are stationary, i.e., they do not vary with vary- 
ing temperature or varying VoF, compared to the BS-based SDM. This 
requirement is also not quite met, which is why probably the JBS-based 
SDM shows a slightly higher systematic error. 

If the methods are evaluated for a data set with varying flow D „ırpı, 
as shown by the results in Fig. 6.7, the situation is quite different. While 
the spline-based methods still work well, the PCA-based SDM is no longer 
applicable because the requirement—a constant time-delay difference— 
is not met at all. Contrary to the PCA-based SDM, a visual comparison 
of the spline-based methods with each other is hardly possible in this 
presentation form. Only for high absolute time delays, corresponding 
to low temperatures in the water data sets, where air bubbles are more 
often present in the experimental setup, the JBS-based approach shows 
more outliers, but further evaluations are needed for a general statement. 


138 


6.5 Signal-dynamics method 


PCA-based BS-based JBS-based 
0.10 T T | | T | 
> 0-05 a 4 - 4 
= 0.00 N/INVN 
= _0.08 | y L | 
0. | | L i 


10 | i L | | 
68 69 70 71 72 68 69 70 71 72 68 69 70 71 72 
t/us t/ps t/us 


Figure 6.8 Estimated interfering signals for the data set D,.tpı- The different interfering 
signals result from different signal packages, which were created using a signal dynamic 
of 50 ns with an adaptive overlap as described in the introduction to Section 6.5. 


6.5.2 Interfering signals compensation 


In the alternative approach, the interfering signals are first determined 
and removed, and then the time-delay differences are estimated using 
a state-of-the-art (SotA) method, instead of estimating the time-delay 
differences directly. As already mentioned in Section 4.2.3, the ZCM 
according to Kupnik et al. [63] is used for the subsequent TDE after the 
interfering signals compensation. In Figure 6.8 all estimated interfering 
signals are shown that result when the interfering signals are estimated 
separately for each of the signal packages. Equally to the evaluation 
of the direct TDE, the data set is D,.ırpı, the signal dynamics are set 
to 50 ns per package, the point clouds are concatenated for the spline- 
based approaches, and the PCA-based SDM uses the constant amplitude 
assumption. Although most signals are estimated to be quite similar 
for all methods, some differences occur that cannot be explained by the 
nonstationarity of the interfering signals. The JBS-based approach, which 
is expected to be the most robust method, shows the least variation of 
the estimated interfering signals. Therefore, it is likely that the reasons 
for the variation are the noise influence and outliers in the point clouds 
that deteriorate the individual BS-based shape fit of the point cloud. 
Furthermore, an incorrect assumption of the signal model, e.g., that the 
amplitude of the target signals does not depend on the temperature, 
could be the cause. 
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Figure 6.9 TDE results using the local interfering signals compensation for the undamped 
water data sets with constant and varying VoF. Overlapping TDEs between the signal 
packages are resolved by averaging. 


If each signal package is compensated by the respective estimated 
interfering signals, application of the zero-crossing-based TDE yields 
the time-delay differences for all signals in the package—this is called 
the local interfering signals compensation. Any overlapping estimates 
resulting from the overlap of the signal packages can be resolved by 
averaging. The results for the data sets D,.ırp1, Dıyırpı can be observed 
in Fig. 6.9(a) and Fig. 6.9(b), respectively. In Fig. 6.9(a) it can be seen that 
the PCA-based SDM continues to perform well for the constant flow 
data set. However, the BS-based point cloud processing does not work in 
combination with the local interfering signals compensation (red curve). 
The inconsistency at 69.7 jis is caused by a gap in the measurement signals, 
which results from the automatic outlier removal of signals degraded 
by bubbles in the experiment. Figure 6.9(b) shows that the PCA is still 
not comparable to the spline-based methods for varying flow data sets. 
This is due to the fact that the abrupt changes of the flow lead to non- 
compact point clouds, which cannot be correctly processed by the PCA. 
In contrast, the spline-based methods show that they are suitable for 
data sets with highly varying flow velocity because the curved form of 
the BSs is flexible enough to fit even point clouds with large gaps. 
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Figure 6.10 Automatic outlier detection to remove bad estimations of interfering signals 
from Fig. 6.8. For better visualization, an offset has been added to the results from the 
PCA and the JBS-based method. Left: the detected outliers are shown in dots. On the right, 
the final estimated interfering signals can be seen, resulting from averaging after outlier 
removal. 


In the global interfering signals compensation, for each signal pack- 
age the same interfering signals estimate is used for compensation. This 
can reduce the effect of bad estimates of the interfering signals, which 
result from outliers in some of the point clouds. To this end, the different 
estimates are collected and outlier signals are removed by applying the 
outlier detection algorithm proposed in Section 4.4. The remaining esti- 
mates for the interfering signals are averaged to obtain the globally valid 
estimate—in the following, using this globally valid estimate will be re- 
ferred to as global interfering signals compensation. A visual example of 
this procedure can be observed in Fig. 6.10. On the left side, all estimates 
for the interfering signals are plotted for each point cloud processing 
method—an offset has been added for better visualization. While the 
remaining estimates after the outlier detection are drawn solid, all signals 
classified as outliers are drawn dotted. It is easy to observe that due to 
the adaptive histogram the limit that defines which signals are identified 
as outliers is dependent on the variability of the estimates. In the case of 
lower variability, such as for the JBS-based SDM, even small deviations 
lead to a classification as outlier, while a higher overall variability, such as 
for the PCA- or BS-based SDM, does not classify signals as outliers even 
with larger deviations. The averaged results are then depicted on the 
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Figure 6.11 TDE results of the global interfering signals compensation for the undamped 
water data sets with constant and varying VoF, resulting if all measurement signals are 
compensated by subtracting the same averaged interfering signal. 


right side. It shows that the resulting averages of the BS- and JBS-based 
SDM are quite close, leading to the expectation that the subsequent TDE 
based on this estimate will perform similarly. 

Application of the global interfering signals compensation to the TDE 
for the water data sets D,.ırpı, Dıyırpı yields the estimation errors de- 
picted in Fig. 6.11. Through the usage of a global estimate for the interfer- 
ing signals, the SDM now produces similar results for the constant and 
varying flow data sets. Only slightly increased stochastic components 
at measurements with low time-delay differences can be observed (see 
the small ripples in Fig. 6.11(b)). Moreover, the expected similar behavior 
of the BS- and the JBS-based SDMs is confirmed. An important point 
to notice is that even a small phase shift in the estimated noise signals 
(Fig. 6.10 on the right) can lead to large estimation errors. This is shown 
by the larger errors of the PCA-based SDM in Fig. 6.11, whereas the inter- 
fering signals estimate exhibited by the PCA-based SDM is only slightly 
phase-shifted compared to the BS- and JBS-based SDM. 

In summary, each of the SDMs has its own advantages and disadvan- 
tages. The PCA-based point cloud processing, regardless of whether 
direct TDE, local or global compensation is used, can only be applied 
with sufficient estimation quality to constant flow data sets. For data 
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sets with varying flow velocity, the spline-based methods are generally 
better. In addition, the spline-based SDM using the global interfering 
signals compensation has been shown to have the lowest estimation error, 
regardless of whether the method was applied to a constant or a varying 
flow data set. 


6.5.3 Hyperparameter optimization and overall results 


The previous two sections only represent a proof-of-concept of the SDMs 
with a standard set of hyperparameters using the water data sets as an 
example. While the results could be used to make simple statements 
about the different approaches, they were neither evaluated quantita- 
tively, nor tested on all data sets with all hyperparameter combinations, 
nor compared to the SotA methods. Therefore, in this section, the SDMs 
are firstly tested for all hyperparameter combinations, where a focus lies 
on the available signal dynamics and their interaction with the available 
SNR. Finally, all SDMs in their best hyperparameter configurations are 
applied to the water and ethanol data sets, evaluated in terms of MAPE, 
and then compared to the best SotA method—the GCC. 

To begin with, the different hyperparameter combinations are tested. 
Figure 6.12 shows the performance of all SDMs in the form of a star 
plot where the categorical hyperparameters are varied for the water 
data sets D,.yrpı, Diyjrpi- A list of the investigated hyperparameter 
sets are found in Table 6.4. Note that for the PCA-based SDMs only 
the individual point cloud generation is possible since the geometric 


Table 6.4 Combinations of hyperparameters available in the SDMs. 


hyperparam- signal point cloud assumptions 
eter set dynamics generation 
apo" Ar, = 50ns individual constant flow 
ae Ar, = 50ns individual constant amplitude 
a Ar, = 50ns individual varying amplitude 
ors Ar, =50ns concatenated - 
07° Ar, = 50ns individual - 
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PCA-based SDM 


Figure 6.12 MAPE results of all three major point cloud processing methods in direct, 
local compensation and global compensation form for both data sets D1 .;pp, and Djy)7p1- 
The MAPE is coded as the distance from the origin for D,,;pp; and Dj yp, in green and 
blue, respectively. Each angle in the star graph represents a hyperparameter configuration 
of the respective processing method. The baseline in its best configuration (best SotA 
method) is drawn dashed. For better visualization MAPE values are capped to 7.5 %. 


interfering signals estimation requires the PCA of two distinct point 
clouds to triangulate the interfering signals. Since the computational 
effort to calculate the hyperparameter combinations full factorial is too 
high, the signal dynamics are fixed to 50 ns in this first hyperparameter 
investigation. The subsequent signal dynamics evaluation will show that 
this setting is justified. Furthermore, the BS- and JBS-based SDMs share 
the same hyperparameter sets, denoted by oar because they only differ 
in the way the splines and the misalignments are estimated, but not in 
the way the point clouds are generated. The results from Fig. 6.12 once 
again shows, that the constant flow data set is a lot easier to manage and 
that the PCA cannot be applied in the varying flow case with sufficient 
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Figure 6.13 Interaction between used signal dynamics Ar, and measurement noise for 
the data set D,.ırpı- Two different noise levels were examined: the original SNR (85.3 dB, 
solid lines) and the deteriorated SNR by adding artificial AWGN (62.9 dB, dashed lines). 


accuracy. Nevertheless, the varying amplitude assumption (0; ““) shows 
slightly better results compared to the constant amplitude assumption, 
which is why it will be used in the following investigations for the PCA- 
based SDMs. For both spline-based approaches, the concatenated point 
clouds are preferable on average. However, even if the hyperparameter 
combination aps is used, both spline-based methods are still better than 
the baseline, except for the direct estimation with individual point clouds. 

Based on these hyperparameter configurations, the influence of the 
available signal dynamics is investigated. To this end, the evaluation of 
each method on the data set D,,)pp; is repeated with different signal 
packaging strategies, where the package size was adapted such that the 
signal packages contain the specified signal dynamics. Starting from 
the minimum 20 ns, the signal dynamics were increased in steps of 5ns, 
10 ns, and 20 ns, depending on the required computational effort of the 


145 


6 Experimental Results 


respective method and the expected variation in the results—for higher 
signal dynamics the step size was consequently increased. Figure 6.13 de- 
picts the results of this investigation. Since the graphic examines several 
aspects of the different SDMs, the results are interpreted step by step. 

Firstly, the MAPE representing the systematic error is investigated. As 
can be easily seen, except for the BS-based SDM using the local compensa- 
tion, the systematic errors are not significantly reduced further for signal 
dynamics greater than 50 ns. The JBS-based SDM would even allow lower 
signal dynamics with only minor degradation of the MAPE. To investi- 
gate the relationship to the available SNR, artificial AWGN is added to 
simulate a lower SNR. The MAPE curve for the lower SNR (drawn as 
dashed line) shows that the lower limit of the signal dynamics, below 
which the SDMs can no longer be used with sufficient quality, shifts 
towards higher signal dynamics. Furthermore, the lower robustness of 
direct TDE is revealed, which is concluded from the fact that the MAPE 
varies much more with different signal dynamics and requires larger 
signal dynamics to perform comparably. Summarizing the dependence 
of the systematic error on the available signal dynamics, it can be said 
that the influence is small as long as the signal dynamics are above the 
lower limit, which is SNR dependent. 

Secondly, the STDSPE representing the stochastic error is investigated. 
Except for the direct estimation form of the SDMs, the stochastic error 
is again hardly dependent on the available signal dynamics. Similar to 
the systematic error, the lower limit of the signal dynamics, below which 
the SDMs can no longer be used, increases. However, additionally to the 
influence on the lower limit, a lower SNR leads to an increased offset of 
the stochastic error. The only exception to these statements is the direct 
estimation. Due to the averaging effect, the PCA-based direct estimation 
SDM can reduce the STDSPE with increasing signal dynamics, but as 
already mentioned this is bought by an increased response time of a 
measurement system based on this algorithm. The other spline-based 
SDMs do not show that they can reduce the stochastic error below the 
thresholds approached by the local, and global compensation methods. 
In summary, as long as the signal dynamics are higher than the SNR- 
dependent limit, further increasing the signal dynamics provides no 
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Figure 6.14 MAPE performance of all SDMs averaged over the different TPs. The results 
achieved with the best SotA method is drawn as a dashed line. The error bars indicate the 
respective best and worst results received if the TPs are investigated individually. 


advantage, except for the direct estimation methods, since they are less 
robust to noise. 

Based on the hyperparameter influences presented above, the final 
overall evaluation of the SDMs is done using signal dynamics of 50 ns, the 
constant amplitude assumption for the PCA-based SDM in its interfering 
signals compensation variants, and the individual point cloud generation 
for the spline-based SDMs. Since the systematic error is more significant 
in the scope of this thesis due to its correlation with the level of the 
interfering signals, the results are compared in terms of MAPE in Fig. 6.14. 
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Here, the averaged MAPE is depicted for the water and ethanol data 
sets without damping mat. Averaged MAPE means that the value is 
determined by averaging the MAPEs of the four transducer pairs per 
data set. Therefore, the basis for comparison—the performance of the 
GCC—is also averaged over all four respective transducer pairs. Note that 
the water data set with damping mat is not included in this evaluation 
since the temperature-dependent damping violates the assumption of 
stationary interfering signals. The results show that the PCA-based SDMs 
cannot be successfully applied to data sets with a strongly varying flow. 
In contrast, using the local or global compensation, the JBS-based SDM 
is significantly better than the best SotA for all data sets. Especially, the 
global compensation works well for both constant and varying flow data 
sets. Although the BS-based SDM is also often better than the GCC, 
some combinations, such as the direct estimation for D,,, or the local 
compensation for D, ,, perform worse. The last point worth mentioning 
here is the range of systematic errors between the transducer pair with 
the largest and smallest MAPE. Some SDMs can significantly reduce this 
range, e.g., the JBS-based global compensation reduces the initial MAPE 
range for the ethanol data set D,. from about 11% to about 3%. 

In summary, except for the PCA-based and some variants of the BS- 
based SDM, all SDMs perform significantly better than the respective 
best SotA method in terms of averaged systematic error. Furthermore, the 
interfering signals compensation does not induce an increased stochastic 
error. 


6.6 Feature-based regression method 


Unlike the SDMs, the FRM requires a training data set, but it also allows 
to obtain an estimate for each signal individually without having to sep- 
arate the signals into packages with some signal dynamics. Since the 
training set should contain varying time-delay differences and different 
process conditions, such as the temperature, the training data set for this 
approach is always chosen from the varying flow data sets. Accordingly, 
the corresponding constant flow data sets, whose time-delay differences 
are not included in the training data set, are used as test data sets to 
predict the time-delay differences of measurement signals that were not 
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seen by the model during the training steps. For a compact represen- 
tation, a simple notation is introduced, to indicate which data sets are 
used as training and test data set, respectively. Exemplary, the notation 
Dyırpı/Dıcyrpı, denotes that the first TP of the varying flow data set 
Dıy is used for training and the corresponding TP of the constant flow 
data D,, set is used to test the estimation performance. 

However, due to the large number of options to select the features, 
the available regression models and their training strategies, the FRM 
requires the specification of much more hyperparameters than the SDM, 
where only the signal dynamics and the point cloud generation are impor- 
tant. A full factorial evaluation takes to much computational effort. For 
this reason, starting bottom-up from the regression model, the hyperpa- 
rameters are investigated consecutively and then fixed to their respective 
best choices in the first two subsections. 

Subsequently to the determination which regression model, train- 
ing strategy, and feature selection method to use, the properties and 
performances of the FRM are evaluated in different scenarios. Firstly, 
the robustness against AWGN is tested, followed by an investigation if 
frames—in this case an overcomplete AWPT tree structure—can further 
improve the results. Finally, the overfitting and the generalizability are 
evaluated, since these properties are highly important in learning-based 
approaches. 

All performances of the FRMs in this section are expressed in terms of 
the MAPE (6.7) for the systematic and the STDSPE (6.10) for the stochastic 
error in order to reduce the entire evaluation of a specific hyperparameter 
set to a single value. 


6.6.1 Selection of the best regression model 


To begin with, the different regression models are investigated. However, 
since the evaluation of the FRM depends not only on the regression model 
but also on the feature selection, the regression models are evaluated 
according to how well they perform on average when using different 
feature selections. Thereby, not all feature selections can be considered 
in the average evaluation, because there are too many combinations of 
feature selections possible, i.e., the quality measure (F-test, PCC, ...), the 
SS or the IFS, and the number of remaining features or thresholds can 
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Table 6.5 Average MAPEs using the three different regression models and their corre- 
sponding training strategies. The averaged MAPE values are obtained by using 
different numbers of selected features from K, = 10... 100. Here, the data sets 
Diyrp1/Picjrpi are used for the training and the test. 


(a) Supervised. 
E{MAPE} Training strategies 
eu) LS WLS WM a 
in % kriging 
Model (5.35a) | 0.3(0.02) 1.5(0.35) 1.3 (0.39) - 
Model (5.35b) | 0.3 (0.04) 1.3(0.56) 3.8 (0.86) = 
GPR (5.35c) - - - 4.7 (3.53) 
(b) Unsupervised. 
E{MAPE} Training strategies 
(MATEN LS WLS WM ap 
in % kriging 
Model (5.35a) | 2.0 (0.14) 1.2(0.49) 1.0 (0.48) = 
Model (5.35b) | 2.0 (0.14) 1.4(0.63) 2.5 (0.42) - 
GPR (5.35c) - - - 2.0 (0.18) 


be combined almost arbitrarily. Therefore, the type of feature selection 
is fixed to the ISS where the parameter of the remaining features K, is 
varied. 

According to sections 5.4 and 5.5, a complete hyperparameter set for 
the FRM contains the following parameters: 


= the feature selection method in combination with the decision 
whether to use the SS or IFS, and a threshold y or top-N value K,; 


= the regression model with its training strategy; 


= the decision whether the regression is to be trained using the 
ground-truth values for the time-delay differences (called the su- 
pervised FRM), or using the estimated time-delay differences from 
the optimization approach described in Sec. 5.5.3 (called the unsu- 
pervised FRM). 
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Additionally, in order to apply the method, the data sets for training and 
testing also need to be specified. 

Table 6.5 shows the resulting average MAPEs of all available regression 
models and training strategies, if the FRM is applied multiple times with 
a differing number of selected features K, and the individual MAPEs 
are averaged. Note that the parameters of the GPR cannot be trained 
using the least-squares or median-based approaches, so simple kriging 
is used here to fit the GPR. To simultaneously evaluate the dependence 
of the results on the number of features K,, the standard deviation over 
the results of the different feature numbers is given in round brackets. 
Comparing the GPR with the linear regression approaches, it performs 
worst in the supervised case and only average in the unsupervised case. 
Therefore, and due to the fact that it requires the highest computational 
effort, this regression approach is not pursued further. Furthermore, it 
is noticeable that while the LS training yields the best results for the 
supervised FRM, the performance in the unsupervised FRM is signifi- 
cantly worse, suggesting that this type of training is less robust to errors 
in the labels. In contrast, the weighted training methods (WLS, WM) 
show similar behavior for both the true labels and the estimated labels. 
A comparison between the regression models indicates that the decision 
whether the model should include an offset or not is only relevant for the 
WM training. The application of this strategy leads on the one hand to 
the best results on average, if combined with the linear model with offset 
(5.35a), and on the other hand to nearly the worst results on average, if 
combined with the linear model without offset (5.35b). Observing the 
standard deviation shows that the dependence on the number of selected 
features is similar for all weighted training strategies, independently of 
the regression model. Only the LS method is significantly less dependent 
on K 

In summary, the use of Gaussian processes for regression is not appro- 
priate for the regression of ultrasonic TDE in this algorithmic context. 
Since the WM method in combination with model (5.35a) performs best 
on average and the LS method has both the least dependence on the num- 
ber of selected features K, and the best results for the supervised FRM, 
these two strategies are further investigated. Here, the LS method is used 
without offset, because the model without offset physically makes more 
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Table 6.6 Evaluation of different methods to select the subset of features to be used for 
the supervised regression. The averaged MAPE values are obtained by using 
different numbers of selected features from K, = 10... 180. The regression 
model (5.35a) was chosen and trained using the WM strategy. As training set 
and test set, D,,;pp; and D,.ırpı were used, respectively. 


Selection method Average MAPE (SS) Average MAPE (IFS) 


ISS 1.3% - 
PCC 0.7% 0.4% 
NCA 1.7% 3.3% 
F-Test score 4.8% 4.6% 
RReliefF 12.8% 17.7% 
Laplacian scores 1.9% 116.4% 
Wrapper 3.4% 15.9% 
LASSO 1.9% 2.2% 
GPR 1.4% 1.4% 


sense when the features are designed to yield zero without time-delay 
difference. 


6.6.2 Feature selection evaluation 


After setting the regression model to the linear model trained by the 
WM or LS method, the different feature selection methods can be quan- 
titatively evaluated. Given the stable results for both supervised and 
unsupervised FRM, all available feature selection methods are tested 
by keeping the regression fixed to the linear model (5.35a) with the 
WM-trained parameters. Again, the different FRM setups are applied 
exemplarily to the water data sets D, ‚ırpı/Dıcjrpı and the results are 
averaged over different top-N values K, to become invariant to this hy- 
perparameter. 

The results of this evaluation are listed in Tab. 6.6. For all methods, 
except for the ISS, where only the subspaces can be selected and not 
the feature individually, the MAPE is given once for the SS and once 
for the IFS. Among the tested methods, only the ISS, the NCA, and the 
Laplacian scores allow feature selection without knowledge of the labels. 
Therefore, these methods are preferable if the MAPEs are comparable. 
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Figure 6.15 Dependence of the best five feature selection methods on the number of 
selected features K,, if the WM parametrized linear regression with offset is applied. The 
regression is trained with the ground truth for the labels. 


Since the RReliefF, the wrapper method, and the Laplacian and F-test 
scores perform significantly worse on average, only the ISS, the PCC, 
the NCA, the LASSO, and the GPR are further investigated. An inter- 
esting fact is that, although the Gaussian processes were not suitable 
as a regression model, the learned length scales in the kernel function 
could be used well as a quality measure for the feature selection. The last 
point that can be seen from this study is that—for the feature selection 
in this specific learning-based TDE approach—the simple PCC is the 
best supervised method and the ISS proposed in this thesis is the best 
unsupervised method. 

In order to determine the number of selected features that yields the 
best overall results, the performances of the five best feature selection 
methods with respect to MAPE and STDSPE are plotted against K, in 
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Fig. 6.15. While the IFS using the NCA or the LASSO becomes unstable 
at high numbers of features, the SS is stable for all methods and all k,. 
However, even though the IFS is unstable for some methods, there are 
also combinations, such as the PCC or the LASSO for a low number 
of features, that result in the best MAPEs that can be achieved overall. 
Since the feature selection should yield a good robustness against noise, 
a low systematic error, and as few selected features as possible to save 
computational effort, the number of features will be fixed to K, = 100 
in the following studies. Figure 6.15 shows that using 100 features is a 
good compromise between the objective to have few features and high 
robustness against noise (see the subfigure at the bottom on the left 
side). The LASSO and the NCA selection methods are not considered 
further, because in some configurations their systematic error is too high. 
Furthermore, the GPR is also discarded due to the significantly larger 
computational effort with only comparable results. 


6.6.3 Robustness against additive white Gaussian noise 


In the previous section, the many possible combinations of feature selec- 
tion methods, regression models, training strategies and label choices 
for the training were evaluated and reduced to a few remaining possi- 
ble combinations. The FRM with one of these remaining configurations 
can now be investigated regarding robustness against noise. Since the 
unsupervised FRM is most comparable to the SotA, which also assumes 
no known labels, the features are selected by the ISS (K, = 100) and 
the linear regression model trained with estimated labels (and the WM 
strategy) is applied to the TDE. Although the present SNR cannot be 
further increased, the performance of the FRM can be investigated for 
decreased SNR by adding artificial AWGN to the training data set, the 
test data set, or both. 

To this end, the water data sets D, ‚ırpı/Dıcyrpı are selected and de- 
teriorated by adding noise. If the level of added AWGN is varied and 
the results of the FRM are plotted over the respective noise level, the 
diagram in Fig. 6.16 is obtained. Depending on which data set the AWGN 
is added to or if the data is smoothed by a prefilter, the different curves 
result. Note that the MAPE of the baseline (AMDF) is not drawn, since its 
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Figure 6.16 Impact of AWGN in the training and test data on the regression accuracy. 
Left: Influence on the systematic error evaluated by the MAPE. Right: Influence on the 
stochastic error evaluated by the STDSPE. The level of the added noise is specified by ø, 
where the subscripts (-)-p, and (-)rs denote whether the noise is added to the training data 
set, the test data set, or both. For reference, the STDSPE of the best SotA method in terms 
of robustness against AWGN (the AMDF) is given. 


systematic error is higher than the visible range. Some interesting effects 
can be observed in Fig. 6.16: 


= Firstly, if the noise is only added during training a slight improve- 
ment of the stochastic error on the test set occurs (see STDSPE 
of the red curve). This can be explained by the fact that this pre- 
sumably shifts the feature selection toward subspaces with higher 
SNR during training. At the same time, however, the systematic 
error increases significantly, which is why the noise level should 
nevertheless be kept as low as possible in the training set. 


= Secondly, a bad SNR in the training data set is worse for the sys- 
tematic error than a bad SNR in the test data set, which can be 
concluded by the MAPE of the red and green curve compared 
to the MAPE of the blue curve, where the noise is only added 
to the test data set. The purple curve shows that prefiltering of 
the measurement signals using a digital filter can alleviate this 
problem, but this is only effective up to a certain limit. 


= Lastly, even without prefiltering, the stochastic error is always 
smaller than the best SotA method. Note that although not drawn 
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here, the systematic error of the AMDF method is almost indepen- 
dent of the added noise level, leading to the fact that the FRM is 
only better than the SotA if the SNR during training is sufficient. 


6.6.4 Application of frames in the feature generation 


As the evaluation of the different feature selection methods and the de- 
pendence on the hyperparameter K, has shown, the features, which are 
available for the regression, are an important influencing factor. There- 
fore, not only the selection process is relevant, but also the initial set of 
all features. Remembering Section 5.3, the features used in the regres- 
sion are created by the absolute and phase differences of the coefficients 
that result from the AWPT of the measurement signals. However, the 
AWPT also requires the filter types and tree structure to be uniquely 
defined. Since the two preliminary studies [151, 159] have shown that 
the 6th level of the AWPT tree provides a good compromise between 
time and frequency resolution, all previous FRM evaluations are based 
on the coefficients of the 6th level in the tree. Nevertheless, the valida- 
tion of this choice is verified again in this section. To this end, the set 
of available initial features is extended by using all coefficients of the 
5th, 6th and 7th level in the AWPT to generate the features. This has 
the special effect that the subspaces covered by the coefficients are now 
overlapping in the time-frequency domain, which induces a redundancy 
in the representation—also called frame. Furthermore, the information 
is redundant, since, e.g., the coefficients from the 7th level contain all 
the information that is already contained in the coefficients from the 6th 
level. 

The bar graph in Figure 6.17 visualizes the MAPEs—averaged over 
all TPs and the respective best and worst results—for the differently 
configured FRMs both when using only the 6th level coefficients and 
when using the redundant coefficients of the 5th, 6th and 7th level, which 
build a frame. In order to make the set of selected features comparable, 
the number of features is set to K, = 100 for all selection methods. 
Except for the ISS in combination with the LS regression and estimated 
labels, using frames yields similar or worse results. There are several 
possible explanations for this behavior, which have to take into account 
the different characteristics of the feature selection methods, the choice 


156 


6.6 Feature-based regression method 


WM regression with offset LS regression without offset 
T T T T j : 
3 | 3H Es, = 6 [| 
E <, < {5,6,7} 
= 
g 2} J 2b 4 
A 
< 
S: gll a aN | 
0 | | i | 0 | T a | 
S g S & 
> © rn >To Sd 
oF ELAS NS PS OL NF ONL 
& gz Ka Y e £ £ Ê Q E 
S L Ss Ss S S S F 
N N 


Figure 6.17 Average MAPEs over all TPs of the data sets D,,,/D,, when using the 6th 
level AWPT coefficients (s,, = 6) to generate the features, and when using all coefficients of 
the 5th, 6th and 7th level (s,, € {5, 6, 7}) to generate the features. The error bars indicate 
the respective best and worst result received if the TPs are investigated individually. For all 
feature selection methods, K, = 100 applies. 


of the labels, and the training strategy of the regression. Firstly, the PCC 
is considered, where the features are selected according to how well 
they are linearly correlated with the ground truth. Here, the results 
are very similar regardless of regression training or selection strategy, 
which in this context means, whether to use SS or IFS. This suggests that 
the selected features of the 6th level already contain the best features, 
confirming the validity of the initial choice of coefficients for the feature 
generation. In fact, a detailed examination of the resulting selection 
masks (not shown here due to the poor visual evaluability) reveals that 32 
features between the initial feature set with s,, = 6 and with s,, € {5,6,7} 
are equal. 

Next, the ISS is considered, where after several iterative selection steps, 
the features with the largest variances are selected. Here, the application 
of frames for the feature generation leads to worse results when the WM 
regression is used, and to better or similar results when the LS regression 
is used (cf. the red and blue bars with the ISS labels in Fig. 6.17). Two 
interesting conclusions can be drawn from this behavior. Firstly, even 
though the LS regression is not very robust to errors in the labels, as 
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shown in Sec. 6.6.1, it can better balance different feature selections by 
suppressing features with low or poor influence in regression. This state- 
ment is based on the fact that the MAPEs of the LS-based regression 
with ISS and supervised training are significantly better than the MAPEs 
of the WM-based regression. Secondly, the ISS returns a worse perform- 
ing selection of features when the initial feature set is larger due to the 
frames. Obviously, the objective to search for the maximum variance 
cannot perfectly represent the quality measure that as much information 
as possible should be contained in the selected subspaces. 

The conclusion of this investigation is that in the problem at hand 
the use of frames does not improve the performance, since a sufficient 
number of good features can already be derived from the coefficients of 
the 6th level given the sampling rate 50 MHz. 


6.6.5 Generalizability of the learning approach 


Since it is the nature of learning-based methods that the quality of the 
training data set significantly affects the performance of the prediction, 
this section presents a quantitative assessment of the generalizability 
of the FRM. The effect that the systematic error on the training data 
set is much lower than that on the test data set is called overfitting. It 
indicates that the model may be too complex and, therefore, can perfectly 
interpolate the training data but hardly predict values not contained in 
the training data. The opposite of this, which means that the performance 
on the test data set is comparable to the performance on the training data 
set, is called generalizability. 

For this reason, first, the performance on the test data set is compared 
with the performance on the training data set, and second, the resulting 
MAPEs are examined when the training data is restricted to include only 
a subset of the process conditions. 


Performance on training data set vs. test data set 


Figure 6.18 shows both the MAPE of the approximation of the trained 
FRM on the training set and the prediction error on the test set. In all 
selection and training scenarios (unsupervised/supervised), the perfor- 
mance on the training data set is worse than or similar to the performance 
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Figure 6.18 Average MAPEs over all TPs with the WM-trained linear regression model 
(5.35a). Blue bars: the MAPEs calculated on the training data set D,,,. Red bars: the MAPEs 
calculated on the test data set D,,. From left to right, the hyperparameter combinations for 
feature selection and the used labels are plotted. The error bars indicate the respective best 
and worst results received if the TPs are investigated individually. For all feature selection 
methods, kK, = 100 applies. 


on the test data set. While the MAPEs are quite similar for the PCC-based 
selection, the ISS yields a smaller systematic error on the test sets. From 
these results, two conclusions can be drawn. Firstly, the linear regression 
model is probably too simple to accurately approximate the time-delay 
differences on the training data set. This can be explained by the fact that 
the time-delay differences of the training data set are highly varying and 
have three clustering points due to the design of the experiment. Both 
properties make the training data set more complex than the test data 
set. Secondly, the PCC-based selections yield in general better feature 
sets, which can also be applied to the estimation of time delays that are 
closer to the boundaries of the process conditions that the FRM was 
trained with. However, this is to be expected since the PCC directly ranks 
the features according to their correlation with the labels and the ISS 
can only make use of quality measures that work without knowledge 
of the labels, such as the signal energy or the variance. Whether the 
ground truth is used for the training of the regression or the labels are 
estimated using the optimization approach seems to have no influence 
in this investigation. 
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Generalizability to different process conditions 


In contrast to the investigation of the MAPE difference between the train- 
ing and test data sets, the influence of constraining the training sets is 
investigated in Fig. 6.19. The aim is to evaluate the extent to which the 
FRM can be used for TDE if the process conditions of the test set are out- 
side the process conditions used for training. To this end, the MAPEs are 
given for the subsets of the test data set, which only contain process condi- 
tions that are also included in the training data set, and for the entire test 
data set. In the following, the MAPEs on the subsets of the test set, which 
contain only the process conditions present in the training data set, are 
referred to as local MAPEs whereas the MAPEs on the entire test set are 
referred to as global MAPEs. Note that, if range{t,}/range{t, } = 1, the 
process conditions of the training and test data set are identical, resulting 
in the local and global MAPEs yielding the same values. Furthermore, 
the results show that the LS-trained regression is less generalizable to 
other process conditions. This can be inferred from the global MAPEs, 
which demonstrate degradation in performance relative to the respective 
best SotA method if the training data set is constrained to contain less 
than half of the available absolute time delay variations. In contrast, the 
local MAPEs are equally good for all regression methods. Note that the 
scale of the plots for the LS-based regressions goes up to 9% due to the 
larger systematic errors. An interesting fact is that the unsupervised 
WM-trained regression shows little dependence on the available train- 
ing set. This can be seen from the behavior of the average and standard 
deviation of the MAPEs, which show a similar magnitude regardless of 
the available training data set. 

Insummary, the supervised WM-trained regression generalizes best to 
different process conditions. Furthermore, even though the unsupervised 
WM-trained regression has a higher variance, it can be applied with 
only little degradation if the available training set contains less varying 
process conditions. However, note that WM-trained regression is better 
than the best SotA method—here the GCC—+regardless of the available 
training set or whether the labels used for training are the ground truth 
or estimated by the approximation approach. 
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Figure 6.19 Performance of the FRM on the test data set D,.ırpı (blue), if the range of 
T, that is available during training is constrained. In red the performance on the local test 
data set Dia = {D C Dycipp1| Tam € t+}, which only contains process conditions that 
are contained in the constrained training data set, is given. In order to get the constrained 
training data sets, different subsets are sampled from the entire training data set. On the 
x-axis, the ranges of the absolute time delays 7, m, which are contained in the used subset 
of the training data set, normalized by the maximum range of the entire data set Dj y7p1 
are given. Since not only the range of the available 7, m but also their location influences 
the result, subsets with equal range are sampled at ten different locations each. The result 
of each used subset is represented by a blue cross mark for the performance on the entire 
test data set and a red dot mark for the performance on the local test data set. The solid 
lines represent the respective average values and the shaded area represents the one-sigma 
range, both calculated based on the results of the ten subsets per given range. The LS- and 
the WM-trained regression are combined with the linear model without offset and with 
offset, respectively. For each method, the features are selected by the ISS with K, = 100. 
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6.6.6 Domain adaption to different media and attenuation 
characteristics 


In the previous section, the generalizability to other process conditions 
was discussed. However, another interesting question is whether the 
training using the data of a specific pipe and medium can also be trans- 
ferred to other media or to other attenuation characteristics. For example, 
the medium can be changed to ethanol or a damping mat can be installed. 
In the field of machine learning, this process of transferring the trained 
model to other domains—different media or attenuation situations in 
this thesis—is called domain adaption. Usually, the trained parameters 
are adapted or fine-tuned to the new domain by considering a small set 
of labeled measurement signals of the target domain. Since this would 
require further calibration measurements, which are not desirable, this 
section investigates the extent to which an already trained FRM can be 
applied to other media and damping properties without fine-tuning or 
parameter adaptation. 

To this end, Figure 6.20 shows the MAPEs of the unsupervised FRM 
that result when different combinations of the available training and test 
data sets are used. Note that for the training data sets, only the varying- 
flow data sets, and for the test data sets only the constant-flow data sets 
are used, since varying time-delay differences have to be present in the 
training data sets. Due to there being three varying-flow data sets and 
three constant-flow data sets, nine different combinations of training and 
test data sets exist. Because the previous investigation only considered 
the water data sets without a damping mat, the necessary initial feature 
set may deviate. Therefore, the results are given for different AWPT tree 
levels. 

It is easy to realize that a transfer from water to ethanol and vice versa 
is not possible (see the bottom row and the right column of Fig. 6.20). In 
contrast to that, regardless of whether the FRM is trained on the water 
data set with or without a damping mat, it can be applied to all water 
data sets with similar quality. The explanation for this behavior is that the 
different SoS means that the learned time ranges cannot be transferred to 
other media. The fact that the unsupervised FRM still shows significantly 
lower MAPEs for ethanol than the best SotA method (10.8 %) shows that 
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Figure 6.20 Domain adaptability of the unsupervised FRM, using the ISS (K, = 100) and 
the WM-trained linear regression model (5.35a). Here, the MAPE values of the different 
TPs using the initial feature sets of the 4th to 9th AWPT tree levels are investigated. In order 
to test the transferability of the training to situations with different media or attenuation 
characteristics, all available training data sets are combined with all available test data sets. 
After training, the FRM is not fine-tuned to adapt to the other domain, i.e., it is evaluated 
to which extent the parameters trained on a specific data set can be applied for the TDE 
of a data set with a different medium or attenuation characteristic. For the baseline, the 
MAPE values, when test and training are performed on the same medium and with the 
same attenuation characteristic, are also included and shown on the diagonal. The data set 
identifiers can be found in Table 6.1. 
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the proposed method also works for other media—provided it is also 
trained on the other medium. 

Another interesting parameter to observe is the AWPT-tree level used 
to generate the initial feature set. For the present ultrasonic system with 
its sampling rate and bandwidth, the AWPT-tree level does not have 
too much influence on the results as long as it is chosen smaller than 
9. If the level is chosen too high, the resulting basis function will be 
too extended in the time domain, which will prevent the resolution of 
different propagation paths. 

The last point to notice is the MAPE of the combination D,,/D;.., 1.e., 
using water as medium ina pipe with a damping mat installed. Com- 
parison with the SotA shows that the FRM cannot further improve the 
systematic error. One possible explanation is the complexity of this data 
set. Even though the damping mat reduces the interfering signals, the 
temperature-dependent damping leads to a violation of the assumption 
that the interfering signals are stationary. Therefore, a regression with 
constant parameters cannot adequately estimate the delay differences 
over the wide temperature range of the test data set. 


6.7 Estimator performances compared with the 
Cramer-Rao lower bound 


The errors in TDE can be separated into systematic errors and stochas- 
tic errors. While the systematic errors arise due to wrong model as- 
sumptions, interfering signals, or calibration errors, the stochastic errors 
originate from measurement noise, jitter, or in this particular case, flow 
fluctuations since the flow is fully turbulent. The systematic errors, which 
in the considered UFM scenario are mainly caused by the interfering 
signals, have been extensively studied and compared to the SotA in the 
previous sections by determination of the MAPEs. Therefore, this section 
focuses on the stochastic errors of the SDM and FRM. To this end, the 
stochastic errors of the respective best SDM and FRM are compared to 
the CRLB—a theoretical lower bound for the variance of unbiased esti- 
mators. Note that the CRLB is a valid theoretical lower bound only if the 
estimators are unbiased and the signal properties such as PSD and ob- 
served time range are accurately known. Furthermore, it only represents 
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a lower limit for stochastic errors that arise due to measurement noise 
and not due to flow fluctuations. Nevertheless, the CRLB can be used 
as a reference point to evaluate the magnitude of the stochastic errors in 
the UFMs. For comparison, the obtained stochastic errors of two SotA 
methods are also examined. 

As the input of the SDM and the FRM are measurement signals with 
the aim to estimate a time-delay difference, both methods can be consid- 
ered time delay estimators. Although they are not unbiased due to the 
existence of interfering signals, the stochastic error can nevertheless be as- 
sessed by removing the systematic error beforehand. Similarly, the ZCM 
and the CC method also represent time delay estimators that are biased 
due to the interfering signals and whose stochastic errors can be com- 
pared to the CRLB. These two methods are specifically chosen because 
they are representative of a time-based TDE and a correlation-based TDE, 
respectively. 

For the investigation of the stochastic error, artificial AWGN of different 
levels is added to the measurement signals of the data set D,.ırpı- The 
level of the artificial AWGN is gradually increased in 20 steps to reduce 
the initially available SNR stepwise from about 85 dB to about 63 dB. Note 
that the actual SNR slightly differs for each measurement, since only 
the amplification of the AWGN can be specified and therefore the actual 
signal energy of the AWGN is random around its specified mean value. 
After applying each method under test to the data set at each noise level, 
the stochastic errors are calculated via removing the systematic errors 
with the median-based lowpass similar to (6.10). The only difference is 
that the errors are not processed as percentages but as absolute errors 
Elm] = T,, — Tm, leading to the standard deviation of the stochastic 
error o(T — 7). 

The resulting stochastic errors depending on the SNR are depicted in 
Fig. 6.21 on the left side. The time ranges used in the ZCM and CC are 
calculated according to the algorithm described at the end of Sec. 6.3. 
For the FRM, which is trained on D,,,)pp;, the hyperparameters are as 
follows: ISS (K, = 100) for the feature selection, approximated labels 
for the training, WM-trained linear model with offset for the regression. 
Furthermore, the results of the SDM in its JBS-based global compensation 
form are plotted. It is easy to observe that the ZCM and the CC method 
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Figure 6.21 Considerations on the standard deviation of the proposed estimation methods 


compared with two SotA methods. As a reference, the results were plotted against the 
CRLB for the present TDE problem. 


perform equally well at high SNR, but the CC method can still be used 
at low SNR without much degradation. Since the SDM is essentially 
the same as the ZCM only with a preprocessing step that removes the 
estimated interfering signals, the stochastic errors are determined by 
the stochastic components in the estimated interfering signals and the 
original stochastic errors of the ZCM. From the fact that the SDM and the 
ZCM perform equally it can be shown that the interfering signals com- 
pensation does not add any stochastic components to the measurement 
signals. The second proposed method, the FRM, can even significantly 
improve the stochastic error. This can be explained by two properties 
of the approach. On the one hand, the feature selection leads to a time- 
frequency selection mask that considers in sum a larger time range for 
the estimation of the time-delay differences. On the other hand, the com- 
bination of several time ranges also achieves an averaging effect over 
different propagation paths, which also reduces the influence of the flow 
fluctuation. 

Since the different methods use different time ranges of the measure- 
ment signal, the CRLB is also different according to (6.11). Therefore, the 
results are shown again as a function of the CRLB in Fig. 6.21 on the right 
side. This plot underlines that the FRM not only performs best in terms 
of stochastic errors but also has the smallest distance to the CRLB, which 
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can be seen from the fact that for a fixed SNR, the absolute distance to 
the dashed line is the smallest. The other investigated methods have 
essentially the same Cramér-Rao efficiency. 


6.8 Comparison and discussion 


All presented methods for TDE (including the SotA) differ in their re- 
quirements, such as labeled or unlabeled training sets, signal packages 
with sufficient signal dynamics, or simply a reference and delay signal. 
Furthermore, the many hyperparameters of the methods, which all lead 
to varying results, make an objective and direct comparison of the meth- 
ods difficult. Therefore, the methods are compared with respect to their 
requirements: the number of hyperparameters, the necessary preprocess- 
ing steps, the achievable measurement dynamics, the practicality, and 
the results. Finally, it is discussed which method is considered the best 
overall in the example application presented—the UFM scenario. 

Firstly, the requirements are discussed. The SotA methods used as 
baselines in this dissertation usually only require a single measurement 
consisting of a reference and delay signal. In contrast to that, the SDM 
requires a signal package consisting of multiple measurements with 
varying signal dynamics, i.e., the absolute time delay, the time-delay 
difference, or both have to vary within the signal package. Then, the time- 
delay difference can be estimated for each signal within the package. 
Similarly, the FRM requires such a signal package for the training, for 
which even the ground truth for the time-delay differences must be 
known as labels in the case of the supervised FRM. 

These requirements on the available measurement signals primarily 
influence the practicality and the measurement dynamics of a measure- 
ment system based on the respective algorithm. The SotA methods, the 
FRM and the SDMs in global compensation form only need one measure- 
ment for the actual TDE, and thus the measurement dynamics are not 
limited by the algorithm. In contrast, the SDMs in local compensation 
form or in direct estimation form can only provide estimates after all sig- 
nals of a signal package have been recorded. This delays the availability 
of an estimation result on average by half of the time needed to record a 
signal package with sufficient signal dynamics. Since the signal dynamics 


167 


6 Experimental Results 


are purely process-driven, this can sometimes be very long. Although 
the FRM and the SDMs in global compensation form can estimate a time- 
delay difference individually for each measurement, and thus do not 
limit the measurement dynamics, they both require a training data set 
beforehand for training in the case of the FRM and for the globally valid 
estimation of the interfering signals in case of the SDM. Therefore, in 
practice, a calibration procedure would be performed before use. While 
the unsupervised FRM and the SDM can be calibrated without knowl- 
edge of the ground truth, the supervised FRM needs the ground truth, 
which can be provided either by another measuring device or by prepar- 
ing the time-delay differences in a calibration environment. However, 
both calibration methods are costly, which is why the unsupervised FRM 
and the SDMs are preferable in practice. 

The next distinguishing feature is the number of hyperparameters. The 
FRM has the most hyperparameters of all methods, because the AWPT, 
the feature selection, the regression model and its training strategy must 
be defined in order to uniquely define the FRM. Nevertheless, the inves- 
tigations have shown that only certain combinations of hyperparameters 
lead to very good results, and many hyperparameters, such as the num- 
ber of features have little effect on quality if they are chosen within certain 
limits. This can significantly reduce the number of possible hyperparam- 
eter combinations. In contrast, the ZCM has only two parameters, namely 
which zero crossing is to be used and the time range used for the linear 
approximation of the phase signal. Between the ZCM and the FRM, the 
rest of the methods fall in line, with the correlation-like SotA methods 
having more hyperparameters than the ZC method—namely the win- 
dow type and the processor—and fewer than the SDMs. The SDMs can 
also be differentiated even more finely. Here, the JBS method has clearly 
more hyperparameters since the optimization for the spline fit is more 
complex. 

Finally, the quality of the estimates is compared in terms of systematic 
and stochastic errors. For this purpose, it is assumed that the basic re- 
quirements of all methods are met. Then, the supervised FRM leads to 
the smallest systematic error on average, followed by the unsupervised 
FRM and the SDM in its global compensation form. The SotA methods 
all yield approximately similar systematic errors, which are significantly 
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worse than the SDMs and the FRMs on all data sets. However, this is to 
be expected, as they are not specifically designed to suppress or filter 
interfering signals. While the FRM is a separate learning-based approach, 
those SDMs based on the estimation and subtraction of the interfering 
signals need to be combined with the SotA methods to perform the TDE. 
Although the SDMs are combined only with the ZCM in this dissertation, 
the combination with the other correlation-like methods is also possible, 
but since the results are expected to be similar, this study was not per- 
formed. Regarding the stochastic error, the results show that, again, the 
FRM performs best and the SDM neither improves nor worsens the SotA 
methods. 

In summary, the unsupervised FRM and the SDM with global compen- 
sation are considered to be the best methods because the training data 
can be collected during normal operation, the measurement dynamics 
is not limited and the interfering signals are significantly suppressed 
without the stochastic error becoming worse. 

Here, it should be noted briefly that since the successful application of 
a learning-based approach to TDE, such as the FRM, suggests that TDE 
can also be accurately performed using approaches from deep learning, 
various network architectures were tested in a Bachelor thesis [150] to 
determine the extent to which they can be applied for the present problem. 
However, the thesis has shown that due to the lacking training data 
even the deep neural network trained in a supervised fashion is not 
comparable to the unsupervised FRM, unless a lot of time and effort is 
invested into developing an adapted network architecture and a targeted 
augmentation strategy. For these reasons, deep learning approaches were 
not considered in this dissertation. 
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7.1 Summary 


In this thesis, two TDE techniques for ultrasonic measurement systems 
have been presented, which are specifically designed to filter or be robust 
against interfering signals in scenarios where very small time-delay dif- 
ferences have to be estimated with high accuracy. Since these challenges 
arise in transit-time UFM systems, where the objective is to accurately 
measure the VoF, the experimental evaluation of the proposed TDE meth- 
ods is based on a transit-time UFM setup. For the quantitative evaluation, 
the VoFs in the experiments are converted into the time-delay differ- 
ences that are to be estimated from the ultrasonic measurement signals. 
Based on these ground truths of the time-delay differences, the proposed 
methods can be evaluated in terms of systematic and stochastic error. 

Firstly, a signal model has been presented that describes the TDE 
problem and the origin of the interfering signals in transit-time UFMs. 
Both presented approach classes—the SDM and the FRM—are based 
on this signal model, although their principles differ, as the FRM is a 
learning-based method for the interference-invariant estimation of the 
time-delay differences and the SDM is a method to estimate and remove 
the interfering signals before applying a SotA method for TDE. 

The requirements of the approaches are that multiple measurements 
with varying time delays are available, where each measurement signal 
needs to have a small bandwidth. Furthermore, the additive interfering 
signals have to be stationary, i.e., they do not vary between consecutive 
measurements, and symmetric, i.e., they are additively superimposed 
equally on both the reference and delay signals. Evidence that these 
requirements are met has been provided by a system identification, which 
also quantitatively evaluates the extent to which the assumptions of 
stationarity and symmetry hold. It has been shown that the amplitude 
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and the time delay of the interfering signals are slightly temperature- 
dependent, which limits the accuracy of the interfering signals estimation, 
and thus the achievable reduction of the systematic errors. 

For the separation of the interfering signals, the SDM employs the 
principle that the time delays of the target signals vary when multiple 
measurements are considered, while the interfering signals are station- 
ary. Therefore, this class of methods represents a generalization of the 
methods proposed by Roosnek [99], where it was assumed that the time 
delays in the considered measurement package uniformly cover one pe- 
riod of the fundamental oscillation. The proposed SDM can significantly 
mitigate the boundary condition, so that the time delays can be arbitrarily 
distributed within the measurement package and the variation of the 
time delays only needs to be large enough to outweigh the measurement 
noise. Within the scope of this work, different variants of the SDM have 
been presented, with the JBS-based global compensation performing the 
best in the exemplary UFM scenario. The JBS-based SDM in its global 
compensation form reduced the systematic error on all considered data 
sets to less than half compared to the best SotA methods. This leads to the 
conclusion that although the interfering signals could not be completely 
removed due to their minor non-stationarity, they can be significantly 
reduced using the proposed SDM. 

The FRM uses the AWPT to get a complex-valued time-frequency 
representation of the measurement signals, which is subsequently used 
to generate features that correlate with the time delay. After an adapted 
feature selection, the remaining features are used as input for a regression 
to directly estimate the time delay differences. The experimental results 
have shown that the FRM works on different media and is always better 
than the best SotA method both in terms of systematic and stochastic 
error. Especially for the stochastic errors, the method is closest to the 
CRLB. With regard to the dependence on the available training data, 
the FRM has been shown to estimate the time-delay difference very 
well even outside the process conditions with which it was trained. In 
addition, an extension of the FRM has been proposed that also allows an 
unsupervised training, i.e., the ground truth of the labels for the training 
data set does not need to be known. Thus, the learning-based approach 
is directly comparable to the SotA method in terms of practicality, since 
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a measurement system can easily be adapted to store past measurement 
signals. 

In summary, the objective to reduce the influence of interfering sig- 
nals in ultrasonic TDEs using the UFM as the application example was 
achieved. The residual estimation errors of the proposed methods in the 
experimental results can be explained by the deviation of the signals 
from the assumed signal model in case of the SDM and by the constraints 
imposed on the regression model in the case of the FRM. 


7.2 Outlook 


The presented SDM puts a strong focus on the representation of the 
signal dynamics as point clouds and their processing. However, this 
limits the approach to the fact that only variable time delays can be used. 
Using an alternative approach as a Bayesian estimation problem, where 
the dependence of amplitudes and time delays on process conditions 
are learned as functions, the interfering signals estimation can then be 
generalized to determine those interfering signals which have the highest 
probability, given an observed signal package. 

In contrast, the FRM is limited to small time delays due to the feature 
design—where small means that the transit time is small relative to the 
period duration of the ultrasonic transducer’s operating frequency. This 
limitation can be overcome by defining other features that combine mul- 
tiple scales in the AWPT tree to be unambiguous over larger transit times. 
This in turn requires abandoning the restriction to linear regressions, so 
other regression types have to be considered in the extended case. 

Apart from improving the proposed methods, there are three direc- 
tions for further research: the application of the proposed TDE algorithm 
to other applications such as the ultrasonic liquid level metering, the 
implementation of the improved path-specific TDE for flow profile imag- 
ing in UFMs and the extension of the FRM to other machine learning 
concepts. 

While the application to other ultrasonic measurement tasks falls into 
the category of further validation of the methods, and flow profile imag- 
ing is an extension of UFM, from a signal processing perspective, the 
extension of the methodology to modern learning techniques is most 
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interesting. Here, different neural network architectures can be explored 
that can efficiently estimate a time delay on one hand and reduce the 
influence of interfering signals on the other hand. For this purpose, con- 
cepts have to be developed on how to reduce the requirements on the 
amount of training data as well as on how to design suitable augmenta- 
tion strategies. 
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